Method for analytical jacobian computation in molecular modeling

ABSTRACT

A method for obtaining analytic Jacobians used in implicit integration methods in the computations for the dynamics of a physical system. With this method, the Jacobian with at least twice the number of digits of accuracy as a numerical Jacobian can be computed. This also results in the implicit integration method being more efficient because a smaller number of iterations are required to solve the nonlinear stage equations of the equations of motion, as well as the ability to take larger timesteps. This speedup in computation is very useful in molecular modeling.

CROSS-REFERENCES TO RELATED APPLICATIONS

[0001] This application is entitled to the benefit of the priorityfiling dates of Provisional Patent Application No. 60/245,730, filedNov. 2, 2000 ; and in addition, No. 60/245,688, filed Nov. 2, 2000 ; No.60/245,731, filed Nov. 2, 2000 ; and No. 60/245,734, filed Nov. 2, 2000; all of which are hereby incorporated by reference.

BACKGROUND OF THE INVENTION

[0002] The present invention is related to the field of molecularmodeling and, more particularly, to computer-implemented methods for thedynamic modeling and static analysis of large molecules.

[0003] The field of molecular modeling has successfully simulated themotion (molecular dynamics or (MD)) and determined energy minima or reststates (static analysis) of many complex molecular systems by computers.Typical molecular modeling applications have included enzyme-liganddocking, molecular diffusion, reaction pathways, phase transitions, andprotein folding studies. Researchers in the biological sciences and thepharmaceutical, polymer, and chemical industries are beginning to usethese techniques to understand the nature of chemical processes incomplex molecules and to design new drugs and materials accordingly.Naturally, the acceptance of these tools is based on several factors,including the accuracy of the results in representing reality, the sizeand complexity of the molecular systems that can be modeled, and thespeed by which the solutions are obtained. Accuracy of many computationshas been compared to experiment and generally found to be adequatewithin specified bounds. However, the use of these tools in the priorart has required enormous computing power to model molecules ormolecular systems of even modest size to obtain molecular time historiesof sufficient length to be useful.

[0004] There are two sources of computational complexity for molecularmodeling tasks involving time integration:

[0005] 1. The particular molecular model which is used to describe thelocations, velocities and mass properties of the constituent atoms, theinter-atomic forces between them, and the interactions between the atomsand their surrounding environment; and

[0006] 2. The particular numerical method used to advance the modelthrough time. Time is advanced repeatedly by very short intervals,called timesteps, until a final time has been reached.

[0007] Substantial work has been completed in reducing the computationalload for molecular models, such as the reduction of model complexity byconstraining higher order modes with rigid body assumptions, simplifyingthe model with rigid or flexible substructuring, Order(N) dynamics,efficient implicit solvent models, and multipole methods for the forcefield models (see, for example, U.S. Pat. No. 5,424,963 on thecommercial MBO(N)D software package). Co-pending applications, U.S.Appln. No.______ , entitled “METHOD FOR LARGE TIMESTEPS IN MOLECULARMODELING,” and U.S. Appln. No.______ , entitled “METHOD FOR RESIDUALFORM IN MOLECULAR MODELING,” both of which are filed of even date, claimpriority from the previously cited provisional patent applications andwhich are assigned to the present assignee, and which are incorporatedby reference herein, describe further improvements in molecular modelsand numerical methods.

[0008] The primary reason molecular simulations are so slow is thatcurrent numerical methods require very small timesteps, typicallybetween 1 and 10 femtoseconds (10⁻¹⁵ to 10⁻¹⁴ seconds). Each timesteptaken requires the computation of a new state (position and motion foreach atom) of the particular molecular model, and then computation ofthe new set of forces resulting from the new state. For example,molecular dynamics simulations of the complex behavior of largemolecules, such as the folding of a protein, typically need to cover atime span from at least a microsecond up to several seconds or evenminutes. With techniques currently in common use, this results in therequirement to take 10⁹ to 10¹⁶ timesteps in the computer simulation.The per-step computations of the state, and especially the forces, growvery expensive as the problem size increases. Even with the fastestcomputers available today, months, years or even centuries of computertime are required to solve such problems even for systems of modestsize.

[0009] Heretofore, it has been widely believed by molecular dynamiciststhat these small timesteps are an inevitable requirement of the need tomaintain accuracy in the presence of the very high frequencies to befound in vibrations of molecular bonds. For example, see Leach,Molecular Modelling Principles and Applications, 1996, p. 328;Berendsen, in Computational Molecular Dynamics: Challenges, Methods,Ideas Deuflhard et al. (ed.), Springer, 1999, pg. 18; Rapaport, The Artof Molecular Dynamics Simulation, Cambridge, 1995, reprinted withcorrections 1998, p. 57; and U.S. Pat. No. 5,424,963.

[0010] This common-sense belief is incorrect, however. The computerscience sub-discipline of numerical analysis has produced an extensivetheory of numerical integration for problems in which high frequenciesexist but are of little interest. These problems are termed “stiff”problems (see, for example, Hairer and Wanner, Solving OrdinaryDifferential Equations II: Stiff and Differential-Algebraic Problems,2^(nd) ed., Springer, 1996). In these cases, it is the stability of theintegration method, not the required solution accuracy, which limits thetimestep. Integration methods exist which have unconditional stability,meaning that in theory they can take arbitrarily large timesteps. Suchmethods have a mathematical property called “L-stability.” Onlyso-called “implicit” integration methods can be L-stable, but very fewimplicit integration methods actually are L-stable. Previously citedco-pending U.S. Patent Appln. No._______, entitled “METHOD FOR LARGETIMESTEPS IN MOLECULAR MODELING,” covers the use of implicit and inparticular L-stable integration methods.

[0011] The present invention covers another critical aspect of allowingthe implicit and in particular L-stable integrators to take largetimesteps, namely, more accurate methods for computing requiredcomponents of the implicit integration methods called “Jacobians”.

[0012] But the same problem of high frequency molecular vibration for MDsimulations causes problems for the calculation of Jacobians. Forexample, the repulsive van der Waal's forces that are generated as theelectron orbitals of two atoms begin to overlap must be accounted for ina molecule. This quantum mechanical effect is often approximated by theso-called Lennard-Jones potential (Rapaport, op. cit.), which models therepulsive force as being proportional to 1/r¹³, where r is the distancebetween the centers of the atoms. This is an extreme stiff interatomicforce which is characteristic of molecular dynamics (MD) simulations andposes particular difficulty for any numerical integration scheme used toadvance time during the simulation. As a result and as mentionedpreviously, most prior art MD integration schemes have timesteps limitedby the high frequencies generated by these extremely stiff repulsiveforces. And in particular, the stiffness of the atomic forces greatlymagnify the difficulty of forming certain required Jacobian matrices.Such Jacobian matrices are a necessary ingredient of any stable implicitintegration scheme, such as described in the immediately citedco-pending application.

[0013] All general-purpose simulation codes provide routines to computeJacobians numerically as follows. For a given equation to integrate {dotover (y)}=ƒ(y,t), the desired Jacobian is J=δƒ/δy and is numericallycomputed:

[0014] where $J = \frac{\Delta \quad f}{\Delta \quad y}$ whereΔ  f = f _(y = y₂)−f_(y = y₁) Δ  y = y₂ − y₁

[0015] Note that the perturbation Δy has to be selected to give a goodresult and may be difficult to choose, especially for stiff systems. Incontrast, analytic Jacobians are computed by solving directly, or inthis case algorithmically, for the equation of the desired derivatives.

[0016] Linearized models are regularly produced analytically for simplesystems. Such linearization is usually performed around an equilibriumsolution. It is common in such packages as ACSL (Advanced ContinuousSimulation Language), (ACSL Reference Guide, Mitchell Gauthier andAssociates, 1996), and SPICE (a circuit simulation package), (R.Kielkowski, Inside SPICE, McGraw-Hill, 1998) to perform equilibriumanalysis followed by linearization. Such linearization is performednumerically.

[0017] In contrast, the Jacobian of the present invention representslinearization about an instantaneous solution of the differentialequations (non-equilibrium) and is generated analytically. It should benoted that another prior art approach to generating analytic Jacobiansis to use automated differentiation tools such as ADIFOR (AutomaticDifferentiation of Fortran) (C. Bischof, et. al., ADIFOR 2.0 Users'Guide, Argonne National Laboratory, 1998) that can symbolicallydifferentiate arbitrary equations. These tools could be used toimplement this invention in practice. However, the structure of theequations must be exploited properly to minimize the computations, toavoid errors due to round off and term cancellations, and to avoid“equation swell” which could limit the size of problems solved.

[0018] Rather, the present invention allows for the calculation ofanalytic Jacobians for the effective implicit integration, includingL-stable integrators, of the equations of motion of molecular models.

SUMMARY OF THE INVENTION

[0019] The present invention provides for a method of modeling thebehavior of a molecule. The method has the steps of selecting a torsionangle, rigid multibody model for said molecule, the model havingequations of motion; selecting an implicit integrator; and generating ananalytic Jacobian for the implicit integrator to integrate the equationsof motion so as to obtain calculations of the behavior of the molecule.The analytic Jacobian J is derived from an analytic Jacobian of theResidual Form of the equations of motion and is described as:$\begin{matrix}{{J = \quad {\begin{pmatrix}\frac{\partial\overset{.}{q}}{\partial q} & \frac{\partial\overset{.}{q}}{\partial u} \\\frac{\partial\overset{.}{u}}{\partial q} & \frac{\partial\overset{.}{u}}{\partial u}\end{pmatrix}{\underset{}{\bigtriangleup}\begin{pmatrix}J_{qq} & J_{qu} \\J_{uq} & J_{uu}\end{pmatrix}}}};{and}} \\{J_{qq} = \quad {\frac{\partial\overset{.}{q}}{\partial q} = {{\frac{\partial({Wu})}{\partial q}\quad {and}\quad J_{qu}} = {\frac{\partial\overset{.}{q}}{\partial u} = W}}}} \\{J_{uq} = \quad {\frac{\partial\overset{.}{u}}{\partial q} = {{- M^{- 1}}\frac{\partial{\rho_{u}\left( {q,u,z} \right)}}{\partial q}\quad {and}}}} \\{J_{uu} = \quad {\frac{\partial\overset{.}{u}}{\partial u} = {{- M^{- 1}}\frac{\partial{\rho_{u}\left( {q,u,z} \right)}}{\partial u}}}}\end{matrix}$

[0020] where q are the generalized coordinates, u are the generalizedspeeds, W is a joint map matrix and M is the mass matrix and ρ_(u) isthe dynamic residual of the equations of motion, and z is−M⁻¹ρ_(u)(q,u,0). The method can also be used for any physical systemwhich can be modeled by a torsion angle, rigid multibody system.

[0021] The present invention also provides for the computer code forsimulating the behavior of a molecule, or any physical system, which canbe modeled by a torsion angle, rigid multibody system. A module in thecomputer code with an implicit integrator includes the analytic Jacobianas described above.

BRIEF DESCRIPTION OF THE DRAWINGS

[0022]FIG. 1 is a representational block module diagram of the softwaresystem architecture in accordance with the present invention;

[0023]FIG. 2 illustrates the tree structure of the multibody system ofthe molecular model according to the present invention;

[0024]FIG. 3 illustrates the reference configuration of the FIG. 2multibody system;

[0025]FIG. 4A illustrate a sliding joint between two bodies of the FIG.2 multibody system; FIG. 4B illustrate a pin joint between two bodies ofthe FIG. 2 multibody system; FIG. 4C illustrate a ball joint between twobodies of the FIG. 2 multibody system;

[0026]FIG. 5 summarizes general computational steps for the ResidualForm method and Direct Form methods used for the analytic Jacobiancomputation;

[0027]FIG. 6 is a chart which summarizes the general computational stepsfor the analytic Jacobian;

[0028]FIG. 7 is a plot of digits of accuracy versus perturbation to showthe accuracy of analytic Jacobian over the numerical Jacobian.

DESCRIPTION OF THE SPECIFIC EMBODIMENTS

[0029] General Description of the Present Invention

[0030] The numerical method used to advance time in the simulation of amodeled molecular system is called the integrator, or integrationmethod. The integration proceeds by discretizing the governing equationsof motion of the molecular system. In the case of an implicitintegrator, this step results in a set of coupled nonlinear algebraicequations (the “stage” equations) and the solution of these equationsyields the state vector for the molecular system at the next time step.

[0031] The present invention aids the solution of the stage equations.Because the atomic forces vary so rapidly over short distances, it isdifficult to propagate atomic trajectories accurately. Small errors inthe atomic positions lead to gross errors in the satisfaction of thestage equations. Since the Jacobian is used solve the stage equationsiteratively, an inaccurate Jacobian leads to trial solution that are farfrom the correct solution. This leads to retraction of trial solutionsand hinders the simulation.

[0032] Numerical Jacobians may be correct in only half their significantdigits. An analytical Jacobian will often be correct in all but the lastdigit. The benefit of this result is that the integrator can take farfewer time steps to simulate the specified interval, allowing fullexploitation of the theoretical stability of the integration method.

[0033] The ease or difficulty in producing the Jacobian dependscrucially upon the formulation used to produce the governing equations.For instance, global Cartesian formulations produce equations with verylimited explicit coordinate dependence. Producing an analytic Jacobianfor such a formulation is well known. On the other hand, using internalcoordinates (in which each molecular subunit's location is expressedrelative to an earlier subunit's location) as independent variables hasgreat benefits. This method is most valuable for any molecule modeledwith any choice of internal coordinates, and in particular when usedwith protein models or other polymeric molecules using “torsion angles”between the residues as internal coordinates. An algorithm for producingan analytic Jacobian for a system formulated in torsion angles isextremely difficult to devise. However, the present invention achievesthis task.

[0034] The present invention addresses a seemingly intractable problem:production of the analytic Jacobian for a formulation using internalcoordinates, and specifically torsion angles, which is generally thoughtto be impractical. In addition to being more accurate than numericalJacobians, the analytic Jacobians are also cheaper (in computing power)to produce. The present invention also recognizes a key result that theJacobian of the state derivatives can be computed by applying a matrixinverse to the Jacobian of the computed torque method. This resultallows significant savings in computer time and effort to construct thealgorithm.

[0035] Furthermore, this method of producing analytic Jacobians formultibody system formulations using torsion angle, internal coordinateshas not been seen in the general MBS literature. The present inventioncan be used for any torsion angle MBS formulation, which can be appliedto many other disciplines besides molecular simulations, including, butnot necessarily limited to, mechanical systems, robotic systems, vehiclesystems, or any other system which could be described as a set ofhinge-connected rigid bodies.

[0036] Overview of Description

[0037] The preferred embodiment is divided into several sections. Thefirst set of sections describes the MD simulation architecture, themultibody system (MBS) definitions, and Residual Form of the MBSequations for the subsequent descriptions. The next set of sectionsdiscusses the definition of the Jacobian, its role in the implicitintegration method, and the computation of the analytic Jacobian usingthe Residual Form. Also shown is the superior accuracy and performanceof the analytic Jacobian vs. the numerical Jacobian. Furtherefficiencies in the computation of the analytic Jacobian are discussed,specifically, exploiting the rigid body MBS to “contract” the size ofthe Jacobian matrix, and exploiting the topological structure of the MBSto eliminate unnecessary computations.

[0038] To solve ordinary differential equations (ODEs), most of theprior art have used equations expressed in the Direct Form, i.e., {dotover (y)}=ƒ(y,t) (where {dot over (y)} means$\left( {{where}\quad \overset{.}{y}\quad {means}\quad \frac{y}{t}} \right).$

[0039] The equations of motion for a biomolecular system can be castinto this form (and called the Direct Form). In molecular modeling, allprior art known to the present inventors have used the Direct Form. Thatis, {dot over (q)}=Wu, {dot over (u)}=M⁻¹ƒ, where q and u aregeneralized coordinates and speeds respectively, so that conventionalODE solution methods can be applied. However, this requires a matrixinversion of M (representing the mass of the system) at a cost ofOrder(N) to Order(N³) floating point operations (depending on algorithmused, where N is the number of degrees of freedom in the system), sincethe natural form of the equations gives rise to inertial couplingbetween the derivatives of the generalized speeds. That is, theequations are most naturally produced in the form {dot over (q)}−Wu=0,M{dot over (u)}−ƒ=0, where M, the mass matrix, depends explicitly uponthe generalized coordinates q, i.e., M=M(q). This fact requires formingand effectively factoring the mass matrix each time the statederivatives are needed by the integrator in integrating the equations ofmotion over time. The generalized joint map matrix W is block diagonaland, although it is also dependent on the coordinates W=W (q), it doesnot have a significant computational cost.

[0040] In accordance with the present invention, a method for thesolution of the equations of a molecular system is expressed in ResidualForm to bypass the customary step of producing the state derivativesdirectly. The Residual Form method has the following steps:

[0041] 1) Discretization of the solution variables. The specific form ofdiscretization is dictated by the particular implicit integration methodused to advance the molecular model in time. Implicit integrationfollows from the Residual Form. Implicit integration, especiallyL-stable integrators and other highly stable integrators, such asimplicit Euler, Radau5, SDIRK3, SDIRK4, other implicit Runge-Kuttamethods, and DASSL or other implicit multistep methods, also provideother advantages for molecular modeling. See, for example, theabove-cited U.S. Patent Appln. No. ______, entitled “METHOD FOR LARGETIMESTEPS IN MOLECULAR MODELING,” filed of even date. As a particularlysimple example, when used with implicit Euler integration, thediscretization is as follows:

{dot over (q)}=(q _(n) −q _(n−1))/h, {dot over (u)}=(u _(n) −u _(n−1))/h

[0042] where h is the timestep.

[0043] 2) Substitution into the residual equations: $\begin{pmatrix}\rho_{q} \\\rho_{u}\end{pmatrix} = \begin{pmatrix}{\overset{.}{q} - {{W(q)}u}} \\{{{M(q)}\overset{.}{u}} - {f\left( {t,q,u} \right)}}\end{pmatrix}$

[0044] 3) Solution of the resulting nonlinear algebraic equations$\begin{pmatrix}\rho_{q} \\\rho_{u}\end{pmatrix} = 0$

[0045] for q_(n) and u_(n).

[0046] The kinematic residual ρ_(q) compares an estimated {dot over (q)}generated from the implicit integrator to the derivatives computed bythe routines for determining the joints of the molecular model, which isdescribed in greater detail below. The second row of the residual isρ_(u), the dynamic residual, which determines the degree to which anestimated {dot over (u)} satisfies the equations of motion.

[0047] The system mass matrix M and the so-called ‘bias-free hingetorque’ ƒ are both state dependent. The bias-free hinge torque isgenerated by the dynamic residual routine when the calculated {dot over(u)} vector passed to the residual routine is zero. In general, thehinge accelerations are a response to applied forces, joint torques, andmotion-induced effects (such as Coriolis and centrifugal forces.) If thesystem were at rest, and subjected only to joint torques, it would beconsidered in a bias-free state. The real system with its actual inputscan be reduced to a bias-free state by computing mathematically. At thecore of the module 54 is a multibody system submodule 66. The analysismodule 56, which communicates with the physical model module 54 and thevisualization module 58, provides solutions to the computational modelsof the molecular systems defined by the physical model module 54. Theanalysis module 56 consists of a set of integrator submodules 68 whichintegrate the differential equations of the physical model module 54.The integrator submodules 68 advance the molecular system through timeand also provide for static analyses used in determining the minimumenergy configuration of the molecular system. The analysis module 56 andits integrator submodules 68 contain most of the subject matter of thepresent invention and are described in detail below.

[0048] The visualization module 58 receives input information from thebiochem components module 52 and the analysis module 56 to provide theuser with a three-dimensional graphical representation of the molecularsystem and the solutions obtained for the molecular system. Manyvisualization modules are presently available, an example being VMD (A.Dalke, et al., VMD User's Guide, Version 1.5, June 2000, TheoreticalBiophysics Group, University of Illinois, Urbana, Ill.).

[0049] The described software code is run on conventional personalcomputers, such as PCs with Pentium III or Pentium IV microprocessorsmanufactured by Intel Corporation of Santa Clara, Calif. This contrastswith many current efforts in molecular modeling which use supercomputersto perform calculations. Of course, further speed improvements can beobtained by running the described software on faster computers.

[0050] Molecular Model and Multibody System Description

[0051] The integrators described below in the submodule 68 operate upona set of equations which describe the motion of the molecular model interms of a multibody system (MBS). To aid the computation of theintegration methods described in detail below, a torsion angle, rigidbody model is used to describe the subject molecule system, inaccordance with the present invention. Internal coordinates (selectedgeneralized coordinates and speeds) are used to describe the states ofthe molecule.

[0052] The MBS is an abstraction of the atoms and effectively rigidbonds that make up the molecular system being modeled and is selected tosimplify the actual physical system, the molecule in its environment,without losing the features important to the problem being addressed bythe simulation. With respect to the general system architectureillustrated in FIG. 1, the MBS does not include the electrostatic chargeor other energetic interactions between atoms nor the model of thesolvent in which the molecules are immersed. The force fields aremodeled in the submodule 62 and the solvent in the submodule 64 in thebiochem components module 52.

[0053]FIG. 2 illustrates the tree structure of the MBS of a subjectmolecule. The basic abstraction of the MBS is that of one or morecollections of hinge-connected rigid bodies 170. A rigid body is amathematical abstraction of a physical body in which all the particlesmaking up the body have fixed positions relative to each other. Noflexing or other relative motion is allowed. A hinge connection is amathematical abstraction that defines the allowable relative motionbetween two rigid bodies. Examples of these rigid bodies and hingeconnections are described below.

[0054] One or more of the bodies, called base bodies 172, have specialstatus in that their kinematics are referenced directly to a referencepoint on ground 174. The system graph is one or more “trees”. Animportant property of a tree is that the path from any body to any otherbody is unique, i.e., the graph contains no loops. The bodies in thetree are n in number (the base has the label 1). The bodies in the treeare assigned a regular labeling, which means that the body labels neverdecrease on any path from the base body to any leaf body 176. A leafbody is one that is connected to only a single other body. A regularlabeling can be achieved by assigning the label n to one of the leafbodies 178 (there must be at least one). If this body is removed fromthe graph, the tree now has n−1 bodies. The label n−1 is then assignedto one of its leaf bodies 180, and the process is repeated until all thebodies have been labeled. This is also done for any remaining trees inthe system.

[0055] To help maintain the relationship between the bodies, an integerfunction is used to record the inboard body for each body of the system.The inboard body for each base is ground and i, the parent or inboardbody 182 for body k 184, is referred to as i=inb(k). Additionally, thesymbol N refers to the inertial, or ground frame 174. A superscript Orefers to the ground origin (0,0,0).

[0056] The symbol for the vector from one point to another contains thename of the two points. Thus, r^(PQ) is the vector from the point P topoint Q. A vector representing the velocity of a point in a referenceframe contains the name of the point and the reference frame: ^(N)ν^(P).Certain symbols to be introduced later relate two reference frames. Inthis case, the symbol contains the name of two frames. Thus, ^(i)C^(k)is the direction cosine matrix for the orientation of frame k in framei. This symbol refers to the direction cosine matrix for a typical bodyin its parent frame. Thus, ^(i)C^(k)(j) indicates the actual body j inquestion. The left and right superscripts do not change with the bodyindex. This is also true for the other symbols. An asterisk indicatesthe transpose: H*(k), for example. A tilde over a vector indicates a 3by 3 skew-symmetric cross product matrix: {tilde over (ν)}w

νv×w. E _(i) is an i by i identity matrix, and 0 _(i) is a zero vectorof length i and 0 _(i) is an i by i zero matrix.

[0057] Rigid Bodies of the Model

[0058]FIG. 3 illustrates the reference configuration 190 of a sample“tree” of the MBS. More than one tree is allowed. A point of each bodyis designated as Q, its hinge point. For example point Q_(k) 186 is thehinge point for body k 184. A fixed set of coordinate axes isestablished in the inertial frame 198. An arbitrary configuration of theMBS is chosen as its reference configuration 190. While in thisconfiguration the image of the inertial coordinate axes is used toestablish a set of body-fixed axes in each body. In the referenceconfiguration each hinge point Q is coincident with P, a point of itsparent body (or extended body.) For each body, point P is called thebody's inboard hinge point. So the inboard hinge point P_(k) 188 forbody k 184 is a point fixed in its parent body i 182. The inboard hingepoint for each base body is a point O 192 fixed in ground. The expandedview that shown in FIG. 2 more clearly shows that point Q_(k) 186 isfixed in body k 184 and point P_(k) 188 is fixed in parent body i 182.

[0059] The hinge point locations define d(k) 194, a constant vector foreach body, and can also be written r^(Q) ^(_(i)) ^(P) ^(_(k)) . Thevector for body k is fixed in its parent body i. It spans from the hingepoint for body i to the inboard hinge point for body k. The vector d(1)196 spans from the inertial origin to the first base body's inboardhinge point (also a point fixed in ground), and can be written r^(OQ)^(₁) .

[0060] For a body, m(k), p(k), and I _(Q) _(k) (k) define the massproperties of body k for its hinge point Q_(k). These are, respectively,the mass, the first mass moment, and the inertia matrix of the body forits hinge point in the coordinate frame of the body. For a rigid bodymade up of a distribution of particles, the mass properties areconstants that are computed by a preprocessing module. The details ofthese computations can be found in standard references, such as Kane, T.R., Dynamics, 3^(rd) Ed., January 1978, Stanford University, Stanford,Calif.

[0061] Let M(k), the spatial inertia of body k for its hinge pointQ_(k), be given by the symmetric 6 by 6 matrix${M(k)} = \begin{bmatrix}{{\underset{=}{I}}_{Q_{k}}(k)} & {\overset{\sim}{p}(k)} \\{- {\overset{\sim}{p}(k)}} & {{m(k)}{\underset{=}{E}}_{3}}\end{bmatrix}$

[0062] Each joint in the system is described by geometric data. Forinstance, a pin joint is characterized by an axis fixed in the twobodies connected by the joint. The particular data for a joint dependson its type. The number n, the inb function, the system mass properties,the vectors d(k), and the joint geometric data (including joint type)constitute the system parameters.

[0063] Joints and Generalized Coordinates of the Model

[0064] FIGS. 4A-4C illustrate the joint definitions of the preferredembodiment of the MBS: the slider joint 100, the pin joint 102, and theball joint 104. Each joint allows translational or rotationaldisplacement of the hinge point Q_(k) 106 relative to the inboard hingepoint P_(k) 108. These displacements are parameterized by q(k) 110, thegeneralized coordinates for body k. In passing, it should be noted thatgeneralized coordinates are examples of generalized quantities, whichrefer to quantities that have both rotational character andtranslational character. For instance, a generalized force acting at apoint consists of both a force vector and a torque vector. Thegeneralized coordinate q(k) for the slider joint 100 is the slidingdisplacement x 112. The generalized coordinate q(k) for the pin joint102 is the angular displacement θ 114. The generalized coordinate q(k)for the ball joint 104 is the Euler parameters (ε₁,ε₂,ε₃,ε₄) 116.

[0065] Each joint may be a pin, slider, or ball joint; or a combinationof these joints. Many other joint types are possible, including, but notlimited to, free joints, U-joints, cylindrical joints, and bearingjoints. For instance, q(k)=(x, y, z), the inertial measure numbers ofthe vector from the base body inboard hinge point to the base body hingepoint express the base body displacement in ground as three orthogonalslider joints. A free joint consists of three orthogonal slider jointscombined with a ball joint, and has the full 6 degrees of freedom.

[0066] The collection of generalized coordinates for all the bodiescomprises the vector q, the generalized coordinates for the system.

[0067] Given the generalized coordinates for a particular joint, twoquantities: r^(P) ^(_(k)) ^(Q) ^(_(k)) (k), the joint translation vectorand ^(i)C^(k)(k), the direction cosine matrix for body k in its parentare formed. The translation vector r^(P) ^(_(k)) ^(Q) ^(_(k)) (k)expresses the vector from the inboard hinge point P of body k to thehinge point Q of body k, in the coordinate frame of the parent body.Details of these computations depend on the joint type and can be easilyderived. For purposes of this description, access to a function that cangenerate r^(P) ^(_(k)) ^(Q) ^(_(k)) (k) and ^(i)C^(k)(k) given thesystem generalized coordinates is assumed.

[0068] As introduced, the choice of hinge point for each body isarbitrary. However, judicious choice greatly simplifies matters. Forinstance, for pin joints the hinge point should be chosen as a point onthe axis of the joint. For this choice points P and Q remain coincidentfor all values of the joint angle, so the joint translation is zero. Ifthe point Q is chosen at a distance from the axis, points P and Q moverelative to each other:

r ^(P) ^(_(k)) ^(Q) ^(_(k)) (k)=λ×r ^(OQ) ^(_(k)) sin θ−(1−cos θ)( E₃−λλ*)r ^(OQ) ^(_(k))

[0069] where λ is the joint axis unit vector, θ is the joint angle, andr^(OQ) ^(_(k)) is the vector from any point on the axis to point Q.

[0070] For pin joints and ball joints, a point on the axis is alwayschosen as the hinge point. For these joints the translation vector r^(P)^(_(k)) ^(Q) ^(_(k)) (k) is zero.

[0071] For a slider joint, the translation vector r^(P) ^(_(k)) ^(Q)^(_(k)) (k) is q(k)λ.

[0072] The direction cosine matrix for a pin is

^(i) C ^(k)(k)= E ₃cosθ+{tilde over (λ)}sinθ+λλ^(*)(1−cos θ)

[0073] The direction cosine matrix for a slider is E ₃.

[0074] Generalized Speeds of the Model

[0075] Let ^(i)V^(k)(k), the generalized velocity of the hinge point ofbody k measured in its parent i, be parameterized by u(k), a set ofgeneralized speeds. Then: ${{{}_{}^{}{}_{}^{}}(k)} = {\begin{pmatrix}{{{}_{}^{}{}_{}^{}}(k)} \\{{{}_{}^{}{}_{}^{Qk}}(k)}\end{pmatrix} = {{H^{*}(k)}{u(k)}}}$

[0076] Here, the matrix H(k) is called the joint map for this joint. Itis a n_(u)(k) by 6 matrix, where n_(u)(k) is the number of degrees offreedom for the joint (1 for a pin or slider, 3 for a ball, 6 for a freejoint). H(k) can, in general have dependence on coordinates q. Given thegeneralized speeds for the joint, the joint map generates the jointlinear and angular velocity, expressed in the child body frame. Thefollowing are used for the joints: ${{H(k)} = \begin{bmatrix}\underset{\_}{\lambda} & 0 & 0 & 0\end{bmatrix}},{pin}$ ${{H(k)} = \begin{bmatrix}0 & 0 & 0 & \underset{\_}{\lambda}\end{bmatrix}},{slider}$ ${{H(k)} = \begin{bmatrix}{\underset{\_}{\underset{\_}{E}}}_{3} & {\underset{\_}{\underset{\_}{0}}}_{3}\end{bmatrix}},{ball}$ ${{H(k)} = \begin{bmatrix}{\underset{\_}{\underset{\_}{E}}}_{3} & {\underset{\_}{\underset{\_}{0}}}_{3} \\{\underset{\_}{\underset{\_}{0}}}_{3} & {{{}_{}^{}{}_{}^{}}(k)}\end{bmatrix}},{free}$

[0077] The collection of generalized speeds for all the bodies comprisesthe vector u, the generalized coordinates for the system. As before,access to a function that can generate the vector ^(i)V^(k)(k) given(q,u) and a specific joint type, is assumed. Access to a function thatcan compute the derivatives {dot over (q)}(k)={dot over (q)}(q(k),u(k))is also assumed. This routine generates the time derivative of thegeneralized position coordinates:

{dot over (q)}=W(q)u

[0078] where W(q) is a block diagonal matrix that relates {dot over (q)}and u, with each block depending upon the joint type:

[0079] {dot over (q)}=u for pin joint, slider joint $\begin{bmatrix}{\overset{.}{ɛ}}_{1} \\{\overset{.}{ɛ}}_{2} \\{\overset{.}{ɛ}}_{3} \\{\overset{.}{ɛ}}_{4}\end{bmatrix} = {{{\frac{1}{2}\begin{bmatrix}ɛ_{4\quad} & {{- ɛ_{3}}\quad} & ɛ_{2\quad} \\ɛ_{3} & ɛ_{4} & {- ɛ_{1\quad}} \\{- ɛ_{2}} & ɛ_{1} & ɛ_{4} \\{- ɛ_{1\quad}} & {- ɛ_{2\quad}} & ɛ_{4}\end{bmatrix}}\begin{bmatrix}\omega_{1} \\\omega_{2} \\\omega_{3}\end{bmatrix}}\quad {for}\quad {ball}\quad {joint}}$

[0080] where q=[ε₁ ε₂ ε₃ ε₄]* and u=[ω₁ ω₂ ω₃]*

[0081] and a free joint is a combination of 3 slider joints and one balljoint. Note that there are 4 {dot over (q)}'s (derivatives of the Eulerparameters) associated with 3 u 's for ball joints.

[0082] Similarly, ^(i)A^(k)(k), the generalized acceleration of thehinge point of body k in its parent, is given by:${{{}_{}^{}{}_{}^{}}(k)} = {\begin{pmatrix}{{{}_{}^{}{}_{}^{}}(k)} \\{{{}_{}^{}{}_{}^{Qk}}(k)}\end{pmatrix} = {{H^{*}(k)}{\overset{.}{u}(k)}}}$

[0083] It is these generalized coordinates q, and generalized speeds u,the internal coordinates for purposes of this description, of themolecular system which are calculated. Rather than working with thetypical inertial coordinates (x, y, z) and speeds in these inertialcoordinate systems, calculations for the subject molecular system arereduced.

[0084] Calculations of the Equations of Motion

[0085] With the exemplary rigid multibody, torsion angle modeldescribed, the equations of motion can now be calculated. In accordancewith the present invention, the motion of the MBS molecular model isdetermined by the Residual Form. The Residual Form method requirescalculations termed the “first” kinematic calculations to distinguishthem from the “second” kinematic calculations, which are furtherrequired by the Direct Form (which is included in this description forpurposes of comparison).

[0086] First Kinematic Calculations for the Molecular Model

[0087] In the first kinematic calculations, given the internalcoordinates of the molecular system, (q, u, {dot over (u)}) and thesystem parameters, the following position, velocity and accelerationkinematics are computed for each rigid body k of the molecular model.(In passing, it should be noted that when the First Kinematiccalculations are done for the Residual Form method, the {dot over (u)}is passed in as a guess of the solution which the integration methodthen refines to the correct solution. In contrast, {dot over (u)} is setto zero when used for the Direct Form method. This is shown clearly inthe later descriptions of the two methods.)

[0088] For each body k compute:

^(N) C ^(k)(k), r ^(Q) ^(_(i)) ^(Q) ^(_(k)) (k), r ^(OQ) ^(_(k)) (k),^(i)φ^(k)(k),

^(N)ω^(k)(k), ^(N)ν^(Q) ^(_(k)) (k)V(k),

^(N)α^(k)(k), ^(N)α^(Q) ^(_(k)) (k), A(k)

[0089] These computations are done recursively, starting from each basebody and progressing to the leaves.

[0090]^(N)C^(k)(k), the direction cosine matrix for body k in ground isdefined as:

^(N) C ^(k)(1)=^(i) C ^(k)(1)

^(N) C ^(k)(k)=^(N) C ^(k)(i)^(i) C ^(k)(k), k=2, . . . n, i=inb(k)

[0091]^(i)C^(k)(k) comes from the joint routine described above.

[0092] r^(Q) ^(_(i)) ^(Q) ^(_(k)) (k), the position vector from Q_(i),the hinge point of the parent of body k to Q_(k), the hinge point ofbody k, expressed in the parent frame, is defined as:

r ^(Q) ^(_(i)) ^(Q) ^(_(k)) (k)=d(k)+r ^(P) ^(_(k)) ^(Q) ^(_(k)) (k),k=1, . . . n

[0093] r^(P) ^(_(k)) ^(Q) ^(_(k)) (k) comes from the joint routine.

[0094] r^(OQ) ^(_(k)) (k), the position vector from the inertial originO to Q_(k), the hinge point of body k, expressed in the global frame, isdefined

r ^(OQ) ^(_(k)) (1)=r ^(Q) ^(_(i)) ^(Q) ^(_(k)) (1)

r ^(OQ) ^(_(k)) (k)=r ^(OQ) ^(_(k)) (i)+^(N) C ^(k)(i)r ^(Q) ^(_(i))^(Q) ^(_(k)) (k), k=2, . . . n, i=inb(k)

[0095]^(i)φ^(k)(k), the rigid body transformation operator for body k isdefined${{{{}_{\quad}^{}{}_{\quad}^{}}(k)} = \left( \quad \begin{matrix}{{{}_{\quad}^{}{}_{\quad}^{}}(k)} & {{{\overset{\sim}{r}}^{Q_{i}Q_{k}}(k)}^{i}{C^{k}(k)}} \\{\underset{\_}{\underset{\_}{0}}}_{3} & {\quad^{i}{C^{k}(k)}}\end{matrix}\quad \right)},{k = 1},{\ldots \quad n}$

[0096] V(k), the spatial velocity for body k at its hinge point,expressed in the frame of body k, is defined${V(1)}\overset{\Delta}{=}{\begin{pmatrix}{{{}_{}^{}{}_{}^{}}(1)} \\{{{}_{}^{}{}_{}^{Qk}}(1)}\end{pmatrix} = {{{}_{}^{}{}_{}^{}}(1)}}$${{V(k)}\overset{\Delta}{=}{\begin{pmatrix}{{{}_{}^{}{}_{}^{}}(k)} \\{{{}_{}^{}{}_{}^{Qk}}(k)}\end{pmatrix} = {{{{{}_{}^{}{}_{}^{k*}}(k)}{V(i)}} + {{{}_{}^{}{}_{}^{}}(k)}}}},{k = 2},{\ldots \quad n},{i = {i\quad n\quad {b(k)}}}$

[0097] A(k), the spatial acceleration for body k at its hinge point,expressed in the frame of body k, is defined${A(1)}\overset{\Delta}{=}{\begin{pmatrix}{{{}_{}^{}{}_{}^{}}(1)} \\{{{}_{}^{}{}_{}^{Qk}}(1)}\end{pmatrix} = {{{}_{}^{}{}_{}^{}}(1)}}$${{A(k)}\overset{\Delta}{=}{\begin{pmatrix}{{{}_{}^{}{}_{}^{}}(k)} \\{{{}_{}^{}{}_{}^{Qk}}(k)}\end{pmatrix} = {\overset{\_}{A} + {\begin{pmatrix}\overset{\sim}{\omega} & \underset{\underset{\_}{\_}}{0_{3}} \\\underset{\underset{\_}{\_}}{0_{3}} & {2\overset{\sim}{\omega}}\end{pmatrix}\quad {{{}_{}^{}{}_{}^{}}(k)}} + {{{}_{}^{}{}_{}^{}}(k)}}}},{k = 2},{\ldots \quad n},{i = {i\quad n\quad {b(k)}}}$

[0098] where$\overset{\_}{A} = {{{{{}_{}^{}{}_{}^{k*}}(k)}{A(i)}} + \begin{pmatrix}\underset{\underset{\_}{\_}}{0_{3}} \\{{{{}_{}^{}{}_{}^{k*}}(k)}\left( {{{{}_{}^{}{}_{}^{}}(i)} \times {{{}_{}^{}{}_{}^{}}(i)} \times {r^{Q_{l}Q_{k}}(k)}} \right)}\end{pmatrix}}$ ω = ^(k*)(k)(i)

[0099] Of course, the computations can all be computed in a single passif desired.

[0100] After completing these steps for one incremental time step, theMBS can service kinematics requests to compute the (generalized)position, velocity, or acceleration information for any point of anybody. This is done by computing the required information for any pointin terms of the hinge quantities for its body, using standard rigid bodyformulas.

[0101] Residual Computation

[0102] With the first kinematic calculations described above, theresidual computation for the Residual Form method can be determined. Adetailed description of the Residual Form and its application tomolecular modeling is found in the previously cited co-pending U.S.Patent Appln. No._______ , entitled, “METHOD FOR RESIDUAL FORM INMOLECULAR MODELING,” filed of even date. This computation fills in twopartitions of the vector $\begin{pmatrix}\rho_{q} \\\rho_{u}\end{pmatrix}\quad$

[0103] given previously. The first partition is called ρ_(q), thekinematic residual, and the second partition is called ρ_(u), thedynamic residual. The kinematic residual is computed from the differencebetween a {dot over (q)}, which is passed-in from the (implicit)integration submodules 66, and the derivative computed by each joint:

{dot over (q)}−W(q)u=ρ _(q)

[0104] The dynamics residual is also computed. Starting with a givenstate of the molecular model, i.e., given (q,u,{dot over (u)}) and thesystem parameters, a program routine models the ‘environment’ of theMBS. Such routines are readily available to, or can be created by,practitioners in the computer modeling field. The routine takes thevalues (q,u) determined by and passed in from the integration submodules66 and returns (the state-dependent) ${{T(k)} = \begin{pmatrix}{T_{Q_{k}}(k)} \\{F(k)}\end{pmatrix}},$

[0105] the applied spatial force for a body k at its hinge point Q_(k),and σ(k), the hinge torque for the body k. The dynamics residual,ρ_(u)(k), associated with generalized speeds u(k) for the body k is thencomputed by the following steps:

[0106] 1. Perform the calculations for the molecular model by theResidual Form as described above with the passed-in state values(q,u,{dot over (u)});

[0107] 2. Generate {circumflex over (T)}(k), the spatial load balancefor each body k in the model having n bodies:${\hat{T}(k)} = {{{M(k)}{A(k)}} + \begin{pmatrix}{{{{}_{}^{}\left. \omega \right.\sim_{}^{}}(k)}\left( {{{\underset{\underset{\_}{\_}}{I}}_{Q}(k)}{{{}_{}^{}{}_{}^{}}(k)}} \right)} \\{{{{}_{}^{}\left. \omega \right.\sim_{}^{}}(k)}\left( {{{{}_{}^{}{}_{}^{}}(k)} \times {p(k)}} \right)}\end{pmatrix} - {T(k)}}$

[0108] k=1, . . . n

[0109] 3. Compute ρ_(u)(k)

[0110] for k=n to 2 by −1

[0111] ρ_(u)(k)=H(k){circumflex over (T)}(k)−ρ(k)

[0112] i=inb(k)

[0113] {circumflex over (T)}(i)+=^(i)φ^(k)(k){circumflex over (T)}(k)

[0114] end

[0115] ρ_(u)(1)=H(1){circumflex over (T)}(1)

[0116] The Residual Form method evaluates the extent to which the systemdifferential equations are satisfied. Zero residual indicates that theapplied forces are in balance with the inertia forces. However, thisdoes not mean the system is in static equilibrium, but rather that theapplied forces would reproduce the given {dot over (u)} when applied tothe system in the state (q,u). The residuals can be interpreted as thatadditional hinge torque needed to balance the applied and inertiaforces. In the literature this method is known as either inversedynamics, or the method of computed torques. It governs the case wherethe {dot over (u)} are all prescribed. At this point all thecomputations required for the Residual Form are complete. The residualsρ_(q) and ρ_(u) are used directly by the implicit integrator in theintegrator submodule 68.

[0117] Second Kinematics Calculations for the Molecular Model

[0118] To carry out the Direct Form method, calculations in addition tothe first kinematics calculations are required. These additionalcalculations are termed the second kinematics calculations. The valuesP(k),D(k),^(i)ψ^(k)(k),^(i)K^(k)(k) are computed as follows:

[0119] 1. Perform the calculations for the Molecular Model by theResidual Form as described above, i.e., the first kinematicscalculations.

[0120] 2. P(k), the articulated body inertia of each body k, isinitialized.

P(k)=M(k), k=1, . . . ,n

[0121] 3. The objects below are then generated:

[0122] for k=n to 2 by −1

[0123] D(k)=H(k)P(k)H*(k)

[0124] G=P(k)H*(k)D⁻¹(k)

[0125] {overscore (τ)}=E ₆−GH(k)

[0126]^(i)ψ^(k)(k)=^(i)φ^(k)(k){overscore (τ)}

[0127]^(i)K^(k)(k)=^(i)φ^(k)(k)G

[0128] i=inb(k)

[0129] P(i)+=^(i)ψ^(k)(k)P(k)^(i)ψ^(k)*(k)

[0130] end

[0131] D(1)=H(1)P(1)H*(1)

[0132] The functional dependence of these quantities is only upon thegeneralized coordinate q. Therefore, the first kinematics calculationsare programmed in anticipation of performing the second kinematicscalculations.

[0133] Forward Dynamics Calculations

[0134] Finally, {dot over (u)} is calculated by sweeping inboard, thenoutboard, of the molecule:

[0135] z(k)= 0 ₆, k=1, . . . n

[0136] for k=n to 2 by −1

[0137] ε(k)=ρ_(u)(k)−H(k)z(k)

[0138] υ(k)=D⁻¹(k)ε(k)

[0139] i=inb(k)

[0140] z(i)+=^(i)ψ^(k)(k)z(k)+^(i)K^(k)(k)ρ_(u)(k)

[0141] end

[0142] ε(1)=ρ_(u)(1)−H(1)z(1)

[0143] υ(1)=D⁻¹(1)ε(1)

[0144] {dot over (u)}(1)=υ(1)

[0145] δ(1)=H*(1)υ(1)

[0146] for k=2 to n

[0147] i=inb(k)

[0148] δ(k)=^(i)ψ^(k)*(k)δ(i)+H*(k)υ(k)

[0149] {dot over (u)}(k)=υ(k)−^(i)K^(k)*(k)δ(i)

[0150] end

[0151] With the First and Second Kinematics Calculations, and theForward Dynamics Calculations, the Direct Form method is available.

[0152] Direct Form Method for the Equations of Motions

[0153] The Direct Form method takes the current state (q,u) and computesthe derivatives ({dot over (q)},{dot over (u)}) using the abovealgorithms, which are then used by the integration method to advancetime.

[0154] Given: (q,u)

[0155] Compute: ({dot over (q)},{dot over (u)})

[0156] 1. Compute {dot over (q)} using joint specific routine as above

[0157] 2. Perform above First Kinematics Calculations with {dot over(u)}=0

[0158] 3. Generate residuals ρ_(u) as above

[0159] 4. Negate the residuals ρ_(u)=−ρ_(u)

[0160] 5. Perform Second Kinematics Calculations

[0161] 6. Compute {dot over (u)} using Forward Dynamics Calculationsabove

[0162] The Direct Form method produces the hinge accelerations {dot over(u)} in response to the applied forces acting on the system. FIG. 5summarizes the computation steps of the Residual Form method and theDirect Form method.

[0163] Jacobians in Implicit Integration

[0164] The MD equations which model a molecule (such as a protein), areimplemented as a multibody system (MBS). These equations representNewton's Laws and are expressed as a set of differential equations {dotover (y)}=ƒ(y,t). The differential equations are implemented using asuite of Order(N) multibody dynamics methods. To advance the equationsin time, in accordance with the present invention, an implicit method ofnumerical integration is used, in particular, L-stable implicitintegration methods, such as implicit Euler, Radau5, and SDIRK3.

[0165] An important ingredient of this integration process is formationof the Jacobian of the differential equations. This is$J\overset{\Delta}{=}\frac{\partial f}{\partial y}$

[0166] Since the function ƒ is itself computed by an algorithm ratherthan by an explicit formula, the Jacobian computation represents asubstantial amount of work. In the simplest approach, the Jacobian canbe formed numerically by differencing the derivative routine. This is adelicate operation because the quality of the Jacobian is a tradeoffbetween round-off and truncation errors. Typically half the workingprecision in the result is retained by choosing a good perturbation sizein the difference scheme. In practice, though, this is difficult to do.

[0167] However, the structure of the governing equations may beexploited to improve the Jacobian computation. The exemplary multibodydynamics methods illustrate this. The algorithms involved compute exactderivatives, even though numerical methods are used to execute theformulas. The derivatives obtained are in error by amounts that dependupon round off and the conditioning of the multibody system underconsideration. But no approximations are involved at the equation level.

[0168] In general, G, the iteration matrix used in the Newton loop ofthe implicit integrator has the form G=E−αJ, where E is the identitymatrix and α is be some scalar function of the timestep. See thepreviously referenced U.S. Patent Appln. No._______ , entitled “METHODFOR LARGE TIMESTEPS IN MOLECULAR MODELING,” filed of even date, for adescription of implicit integrators. Changes in step size requirerefactoring G , but not reforming J. Reforming J is needed only when theJacobian is needed at a new state. G is used in a linear subproblemwithin a Newton loop. The following is solved:

GΔy=−r(y _(n) ^(i)),

y _(n) ^(i+1) =y _(n) ^(i) +Δy

[0169] where r(y) is the residual function for that particular implicitintegration method.

[0170] As shown later, J has a special structure, which is inherited byG. This means that equation solving with G can be done at reduced cost,if the structure is exploited.

[0171] The quality of the Jacobian affects the ability to solve thenonlinear equations resulting from discretization in the integrator.Failure to solve the Newton loop may require retraction of a trail stepand reduction of the integration time step. The timestep should becontrolled by accuracy, rather than failures in the Newton loop.

[0172] Computation of Analytic Jacobian

[0173] The Jacobian J is a matrix which represents a linearization ofthe equations of motion. Normally, the governing equations for adynamical system are linearized around an equilibrium state, or perhapsa state of steady motion. In this case, the equations are linearizedaround an arbitrary state so all possible contributing terms should bedeveloped. It is customary to describe J in terms of its partitions:${J\left( {q,u} \right)} = {\begin{pmatrix}\frac{\partial\overset{.}{q}}{\partial q} & \frac{\partial\overset{.}{q}}{\partial u} \\\frac{\partial\overset{.}{u}}{\partial q} & \frac{\partial\overset{.}{u}}{\partial u}\end{pmatrix}\overset{\Delta}{=}\begin{pmatrix}J_{q\quad q} & J_{q\quad u} \\J_{u\quad q} & J_{u\quad u}\end{pmatrix}}$

[0174] Structure of J_(qq) AND J_(qu)

[0175] The {dot over (q)} equation is {dot over (q)}=Wu, where thematrix W has a block-diagonal structure. Each block depends upon thejoint type. Pins and slider joints give rise to 1 by 1 identity blocks.A ball joint generates a 4 by 3 block that expresses the Euler parameterderivatives in terms of Euler parameters and body angular velocitymeasure numbers (generalized speeds).

[0176] From the {dot over (q)} equation above, these two partitions ofthe Jacobian matrix are formed:$J_{q\quad q} = {\frac{\partial{W(q)}}{\partial q}u}$ J_(q  u) = W

[0177] These equations are to be interpreted for symbolic purposes only.In practice, there is no need to generate the matrix W explicitly. Thenon-zero block diagonal elements are filled in as discussed in theprevious section on the kinematic residual.

[0178] Structure of J_(uq) AND J_(uu)

[0179] The {dot over (u)} derivatives are more complicated. Since {dotover (u)}=−M⁻¹ρ_(u),$\frac{\partial\left( {{- M^{- 1}}\rho_{u}} \right)}{\partial q}$

[0180] and$\frac{\partial\left( {{- M^{- 1}}\rho_{u}} \right)}{\partial u}$

[0181] must be computed. (Note that ρ_(u) is the residual of thedynamics routine developed earlier). Here a key result from the field ofNumerical Analysis is used to avoid the derivative of the matrixinverse.

[0182] Suppose A(x)y(x)=b(x), then we can write y(x)=A⁻¹(x)b(x). Ify(x₀)=z is known, and the value of$\frac{\partial{y(x)}}{\partial x}_{x = {x0}}$

[0183] must be obtained, we have${\frac{\partial y}{\partial x}_{x = {x0}}} = {{{A^{- 1}({x0})}\frac{\partial}{\partial x}\left( {{b(x)} - {{A(x)}z}} \right)}_{x = {x0}}}$

[0184] where z=A⁻¹b is held fixed when forming the right-hand side ofthe equation. Applied to the described multibody routines,$\frac{\partial\overset{.}{u}}{\partial x} = {{{- \frac{\partial}{\partial x}}\left( {M^{- 1}{\rho_{u}\left( {q,u,0} \right)}} \right)} = {{- M^{- 1}}\frac{\partial}{\partial x}{\rho_{u}\left( {q,u,z} \right)}}}$

[0185] where z

M⁻¹ρ_(u)(q,u,0). This result avoids the computing of the derivative ofM⁻¹, which is a hypermatrix. The matrix inverse is “pulled-out”. In theabove equation, “x” is either the generalized coordinate q or thegeneralized speed u, depending on the Jacobian partition that is to becomputed.

[0186] In summary, to compute the blocks J_(uq) and J_(uu), three stepsare followed:

[0187] 1. Given (q,u) compute z

−M⁻¹ρ_(u)(q,u,0) using the Direct Form method. This simply says tocompute {dot over (u)} from the current state and save it as thevariable z.

[0188] 2. Compute the analytic Jacobian of the dynamics residualroutine. In this step, the matrix$\frac{\partial}{\partial x}{\rho_{u}\left( {q,u,z} \right)}$

[0189] is to be formed. This quantity clearly depends upon the vector zcomputed in Step 1. Note that the numerical value of the residual ρ_(u)in this step is zero for each element since we are computing theJacobian around a consistent solution of the motion equations. Thepartitions J_(uq) and J_(uu) of the Residual Jacobian are obtained bysubstituting “q” and “u” for the “x” above.

[0190] 3. Back-solve the result of Step 2 with the mass-matrix to obtainthe desired block. The back-solve operation is accomplished in theDirect Method routine by processing a residual vector into a {dot over(u)} vector. The Second Kinematics Step only needs to be performed once,since the back-solves are done at the nominal value of the state. Infact, the Second Kinematics routine must have been called in Step 1while computing z, so the variables should still be cached.

[0191] In words, the Jacobian of our derivative routine can be formed byback-solving the Jacobian of our residual routine. The residual Jacobianis much simpler to compute than the derivative Jacobian. Steps 2 and 3above are derived in the following sections.

[0192] The Residual Jacobian

[0193] The computation of the Residual Jacobian is closely related tothe Residual Form method for dynamics, which is summarized here:

[0194] 1. Perform an outboard pass that computes the kinematic data thatdepends upon position and velocity.

[0195] 2. Make a call to the force routine which generates the atomicforces and consolidates them into spatial loads acting on the bodies.

[0196] 3. Perform another kinematic pass that computes accelerationlevel quantities (using passed-in {dot over (u)}), and combines inertiaforces with the spatial loads from Step 2.

[0197] 4. Perform an inward pass that generates the residual at eachjoint. This pass recursively computes the resultant of the (spatial)inertial and applied forces for the ‘nest’ of bodies outboard of thejoint in question. The residual is the projection of the resultant onthe joint's degrees of freedom, given by the joint map data.

[0198] At a high level the residual computation can be considered todepend upon two kinds of forces: ‘motion forces’ and external forces.The motion forces are computed directly by the multibody system. Theexternal forces are available to the multibody system from a forcemodeling routine that computes the various interatomic forces such aselectrostatics and solvents. A similar procedure is followed whencomputing the Jacobian. The multibody system builds the Jacobian of themotion forces, and combine it with the Jacobian of the external forces.

[0199] Partitioning the Force Field into Effects

[0200] There are many forces that may be acting on the molecule, andthese forces may be computed in various intrinsic coordinate frames thatare most convenient for that particular force “effect”. For example,electrostatic terms may be computed using multipole methods andspherical coordinates, covalent terms may be computed in terms oftorsion and bond angles, and solvent forces may be computed in globalCartesian coordinates. During the computation of the Residuals, theseforces are transformed from their intrinsic coordinate frame to the MBScoordinates.

[0201] The same exchange occurs to compute Jacobians. The nativeJacobians in their intrinsic coordinates are be brought into the MBScoordinates. This requires the use of the chain-rule to transformbetween intrinsic and the MBS generalized coordinates. It is importantthat each effect co-computes its function value and Jacobian, becausemany of the same terms are needed for each computation. Each effect istransformed into a set of spatial loads T_(effect)(k), where k is theindex of a generic body in the system. The totality of these effects isgiven the symbol T(k).

[0202] Effect Jacobians Brought into the Residual Jacobian

[0203] At a high level, the residual routine was previously implementedfrom the equation

ρ_(u) =Hφ(MA−T)

[0204] The implementation uses Order(N) methods which are immediatelyobvious from the equation above. In this equation T is a vector ofspatial loads acting on the pivots of the multibody system, where eachelement is a spatial load (a 6-vector composed of one force and onetorque). It actually represents all effects other than inertia loads orpure hinge loads. The term in parentheses represents the load balancefor each body. The first term is the inertia force, the next is thespatial load. M(k)A(k) is the spatial inertia force for a typical body.This is built from the body mass properties and the spatial accelerationof the body pivot. The spatial acceleration is computed before theresidual routine is executed by the Forward dynamics routine. Theoperator Hφ is implemented in a routine that performs an Order(N) inwardpass.

[0205] Even without knowing anything about the details of thecomputation implied by the equation above, the contribution of$\frac{\partial T}{\partial x},$

[0206] the spatial load Jacobian (the effect Jacobian) to$\frac{\partial\rho_{u}}{\partial x},$

[0207] the residual Jacobian, can be immediately inferred:$\frac{\partial\rho_{u}}{\partial x} = {{{- H}\quad \varphi \frac{\partial T}{\partial x}} + \ldots}$

[0208] The . . . refers to terms not involving the effect Jacobian.Again, q or u for “x”, is substituted, depending upon which partition ofthe Jacobian is being computed.

[0209] The role of an effect Jacobian in the residual Jacobian is thesame as the role of the effect in the multibody equations. This meansthat T contributes to the residual ρ_(u) in the same way that a columnof $\frac{\partial T}{\partial x}$

[0210] contributes to $\frac{\partial\rho_{u}}{\partial x}.$

[0211] Both are processed by the same operator Hφ. This is a crucialpoint because it means that no new method is needed for this part of theJacobian computation. (A different method is need to obtain$\left( {A\quad {different}\quad {method}\quad {is}\quad {need}\quad {to}\quad {obtain}\quad \frac{\partial T}{\partial x}} \right).$

[0212] Thus, given the effect Jacobian, its contribution is assembledinto the residual Jacobian by operating on it with the original residualroutine, treating it as a multi-column set of spatial load vectors. Thisis a direct consequence of the linearity of the equations.

[0213] The columns of the residual Jacobian play the same role in thederivative Jacobian routine as the residual vectors play in the ForwardDynamics routine. The dynamics routine performs a back-solve on a datavector it receives, and doesn't need to know what the data is, just whatoperation to perform on it. This applies to all the routines.

[0214] Adding the . . . terms, the chain rule is used to show the wholeequation:$\frac{\partial\rho_{u}}{\partial x} = {{{H\varphi}\left( {\frac{\partial({MA})}{\partial x} - \frac{\partial T}{\partial x}} \right)} + {H\frac{\partial\left( {\varphi \quad z} \right)}{\partial x}} + \frac{\partial({Hy})}{\partial x}}$

[0215] At this level, there are four contributions to the Jacobian: theinertia forces, the spatial forces, and contributions due to changes inthe operators φ and H. The quantities z and y refer to (MA−T), and φz,respectively, which are held fixed while evaluating the last two terms.The numerical values of these terms are already available from theresidual computation. Another observation about the above equation isthat the operator φ depends only upon q, but not u. Thus, this termremains constant while computing the partition J_(uu). Similarly, thespatial load can be split into its constituent effects, some of which donot depend upon u. In general, this means that knowledge of themultibody equations can be exploited to optimize the computations.

[0216] Up to now, what to do with $\frac{\partial T}{\partial x}$

[0217] once the term has been computed has been described, but adescription of how to form the term has not been made. These details arein the following sections.

[0218] Computing the Effect Jacobians and Combining with the ResidualJacobian

[0219] So far a high-level description of the Jacobian computation hasbeen given. It can be seen that the computation has a very algorithmicflavor to it. There are very distinct phases to the task, just as therewere also for the Forward Dynamics routine described previously. There,computation of atomic forces is clearly the bottleneck step. Yet theoverhead in the multibody equations for dealing with forces is fairlysmall. In that case, a call was made to the force routine, and whatoccurs inside the routine was ignored. When it comes to the Jacobian,this aspect is less true.

[0220] A call is still made to obtain the Effect Jacobian, but there isa lot of processing needed before the Effect Jacobian can be assembledinto the Residual Jacobian. The details of dealing with force fields toproduce Jacobians are covered in the next sections and an example ofincorporating electrostatics is developed. All other loads follow asimilar development.

[0221] Electrostatics as an Example Effect

[0222] The basic premise of electrostatics is that the force between twocharged particles is$F_{ij} = {\kappa \frac{q_{i}q_{j}}{r_{ij}^{2}}{\hat{r}}_{ij}}$

[0223] This is a classical inverse square law. F_(ij), the force actingon particle i due to particle j, depends upon the charge of theparticles, attractive for oppositely charged particles, repulsive forlike charges. The symbol {circumflex over (r)}_(ij) is a unit vectordirected from particle i to particles j; r_(ij) is the distance betweenthe particles; k is a unit-dependent constant related to the strength ofthe forces. Of course, the force acting on particle j is equal andopposite to that acting on the particle i. Hence, given a collection ofparticles interacting through Coulomb's Law given above, the net forceacting on each particle is computed by summing the pair-wise forces. Forthis example, the forces are computed in global Cartesian coordinates.

[0224] With the atomic forces in hand, the multibody forces can begenerated. The system of forces acting on the particles of each body isreplaced by a spatial load acting at the pivot of each body. The atomicforces are first expressed in a body-fixed basis, and then shifted tothe pivot using the station coordinates of the particular atom to whichthe force is bound.

[0225] F_(ij) is a vector. The derivative of F_(ij) with respect todr_(ij), small changes in the particle's relative positions is D _(ij),a tensor. It is such that dF_(ij)=D _(ij)·dr_(ij). For Coulomb force thetensor is${\underset{\underset{\_}{\_}}{D}}_{ij} = {\kappa \frac{1}{r_{ij}^{3}}\left( {{\underset{\underset{\_}{\_}}{E}}_{3} - {3{\hat{r}}_{ij}{\hat{r}}_{ij}^{*}}} \right)}$

[0226] Some observations:

[0227] 1. The Force Jacobian is a matrix of size n_(atoms)×n_(atoms).Each element is a 3 by 3 tensor. The (i,j) block gives the derivative offorce on atom i with respect to small changes in the position of atom j.In general, every force model is required to support an intrinsicJacobian method for analytical processing.

[0228] 2. Storage requirements for the Force Jacobian quickly becomeimpractical. This leads to the notion of interface “contraction” wherethe Jacobian of all the forces acting pair-wise between the atoms arereduced or “contracted” to acting pair-wise between the bodies.

[0229] 3. Because the Coulomb force is a pair-wise interaction, eachforce contributes to two blocks in the overall Jacobian. Thus, eachforce is processed at constant cost, and the overall Jacobian iscomputed at a cost proportional to the number of atoms squared, i.e.,Order(N²) This is the same as the computational cost of the forceitself! This is a rather good result for computing analytic Jacobians. Anumerical Jacobian requires a fresh force computation each time anelement of the state is perturbed. This leads to cubic growth, i.e.,Order(N³), in the cost of the numerical Jacobian. Hence, the analyticJacobian is much cheaper to compute as well as more accurate than anumerical Jacobian.

[0230] 4. The computation of the Jacobian is conveniently done in thesame routine that computes the force (co-computation). However, ittypically needs to be done far less often than the force computation.Therefore, a flag can be used to trigger the Jacobian calculation onlywhen needed.

[0231] Coupling to the Displacement Gradient

[0232] Having obtained the (intrinsic) Force Jacobian, it is necessaryto process the Jacobian further. This is due to the fact that themultibody system is formulated using relative coordinates. The chainrule is applied to each atomic force F(k,i) and is called “coupling tothe displacement gradient”. This denotes the global Cartesian force onatom i, which resides on body k.$\frac{\partial{F\left( {k,i} \right)}}{\partial q_{j}} = {\sum\limits_{p = 1}^{nbod}{\sum\limits_{s \in \quad p}{\frac{\partial{F\left( {k,i} \right)}}{\partial{r\left( {p,s} \right)}}\frac{\partial{r\left( {p,s} \right)}}{\partial q_{j}}}}}$

[0233] where r(p,s) is the position of atom “s” on body “p”.

[0234] The first term in the sum selects an element of the ForceJacobian which was just computed. The quantity$\frac{\partial{r\left( {p,s} \right)}}{\partial q_{j}}$

[0235] is an element of the displacement gradient. A typical term givesthe change in an atom's position due to a small change in a generalizedcoordinate. Note that this term is strictly a kinematical quantityhaving nothing to do with the force computation. Thus, the ForceJacobian can be computed once and then continually reprocessed by thechain rule for each coordinate in the multibody system. This steprepresents a matrix vector multiply, since$\frac{\partial r_{j}}{\partial q_{s}}$

[0236] is a column vector with n_(atoms) entries (each a 3 vector), andthe Jacobian is a square matrix n_(atoms)×n_(atoms), where each elementis a 3 by 3 tensor.

[0237] It is possible to improve this computation, since many of theentries in the displacement Jacobian are known to be zero. This is dueto the fact that incrementing a particular hinge does not displace everyatom in the system, but only those outboard of the displaced hinge, andnot disjoint from the hinge in question. For instance, rotating the basebody induces a change in all atoms of the system. But perturbing atorsion angle on any terminal body induces a change only to those atomsresident in the terminal body. Therefore, roughly speaking, about halfthe work may be saved by optimizing the computation. This reductioncomes from a strictly knowledge-based approach.

[0238] Interface Contraction

[0239] In the process of forming T(k), the spatial load on body k, theload comes from the atomic forces acting on atoms that belong to body k.Each force is transformed from global to local coordinates, and thenshifted to the body pivot. A concise statement of this procedure is:${T(k)} = {\sum\limits_{i \in \quad k}{{\varphi \left( {k,i} \right)}{T\left( {k,i} \right)}}}$

[0240] The operator φ(k,i) is${\varphi \left( {k,i} \right)} = {\begin{bmatrix}{\underset{\underset{\_}{\_}}{E}}_{3} & {\overset{\sim}{\rho}\left( {k,i} \right)} \\{\underset{\underset{\_}{\_}}{0}}_{3} & {\underset{\underset{\_}{\_}}{E}}_{3}\end{bmatrix}\begin{bmatrix}{{{}_{}^{}{}_{}^{}}(k)} & {\underset{\underset{\_}{\_}}{0}}_{3} \\{\underset{\underset{\_}{\_}}{0}}_{3} & {{{}_{}^{}{}_{}^{}}(k)}\end{bmatrix}}^{*}$

[0241] where ρ(k,i) are the fixed station coordinates of atom i on bodyk. Note that the new quantity T(k,i) appears. This is just the atomicforce turned into a spatial load at the atomic position of atom i:${T\left( {k,i} \right)} = \begin{bmatrix}{\underset{\underset{\_}{\_}}{0}}_{3} \\{F\left( {k,i} \right)}\end{bmatrix}$

[0242] The first element of the atomic spatial load is zero becausethere is no torque exerted by the force field on individual atoms.

[0243] Now, T(k) relates atomic forces to body spatial loads. So, thederivative of this equation relates differential atomic forces todifferential spatial loads: $\begin{matrix}{\frac{\partial{T(k)}}{\partial q_{j}} = {{\sum\limits_{i \in \quad k}{{\varphi \left( {k,i} \right)}\frac{\partial{T\left( {k,i} \right)}}{\partial q_{j}}}} + {\frac{\partial{\varphi \left( {k,i} \right)}}{\partial q_{j}}{T(k)}}}} \\{= {{T_{1}(k)} + {T_{2}(k)}}}\end{matrix}$

[0244] The second term, T₂(k), in this equation is discussed at the endof this section, as it involves the spatial loads, but not the loadderivatives. This means the term can be treated generically, withoutworrying how the spatial loads were computed.

[0245] Substituting the definition of T(k), into T₁(k): $\begin{matrix}{{T_{1}(k)} = {\sum\limits_{{i\varepsilon}\quad k}^{\quad}{{\varphi \left( {k,i} \right)}{\sum\limits_{p = 1}^{nbod}\sum\limits_{{s\varepsilon}\quad p}}}}} \\{= {\sum\limits_{p = 1}^{nbod}{\sum\limits_{{s\varepsilon}\quad p}{\frac{\partial{T(k)}}{\partial{r\left( {p,s} \right)}}\frac{\partial{r\left( {p,s} \right)}}{\partial q_{j}}}}}}\end{matrix}\frac{\partial{T\left( {k,i} \right)}}{\partial{r\left( {p,s} \right)}}\frac{\partial{r\left( {p,s} \right)}}{\partial q_{j}}$

[0246] where now the symbol$\frac{\partial{T(k)}}{\partial{r\left( {p,s} \right)}}\overset{\Delta}{=}{\sum\limits_{{i\varepsilon}\quad k}{{\varphi \left( {k,i} \right)}\frac{\partial{T\left( {k,i} \right)}}{\partial{r\left( {p,s} \right)}}}}$

[0247] is an element of the row-reduced Force Jacobian. This Jacobianrelates differential (body) spatial loads directly to differentialatomic displacements. Again, r(p,s) refers to the global position of theatom s on the body ρ. The term$\frac{\partial{T(k)}}{\partial{r\left( {p,s} \right)}}$

[0248] is formed by summing each atomic Force Jacobian element into thedestination element in the reduced Jacobian, weighted by the atoms'φ(k,i) matrix. Each element of the row-reduced Jacobian is a 6 by 3matrix. Hence the rows of the Force Jacobian have been contracted. Thecontraction is evident in the notation: the numerator has only a bodyindex, while the denominator has both a body and an atom index.Depending on the number of atoms per body, the row reduction can providea savings in both storage and execution time when differential spatialloads must be formed.

[0249] Note that the row-reduction procedure only needs to be done oncebefore computing the Residual Jacobian. The overhead of performing thereduction is more than offset by the reduced cost of the smaller matrixvector products which must be formed. Note that in forming$\frac{\partial{T(k)}}{\partial{r\left( {p,s} \right)}}$

[0250] there is no need to save the elements of the atomic ForceJacobian. That is, each element$\frac{\partial{T\left( {k,i} \right)}}{\partial{r\left( {p,s} \right)}}$

[0251] only needs to be available while its contribution to$\frac{\partial{T(k)}}{\partial{r\left( {p,s} \right)}}$

[0252] is being computed. So, more than one element of the big ForceJacobian is not required at a time.

[0253] The global coordinates of a typical atom, r(p,s), are computed interms of r(p), the global coordinates of body p's pivot, and ρ(p,s), thestation coordinates of the atom:

r(p,s)=r(p)+^(N) C ^(k)(p)ρ(p,s)

[0254] By differentiating we find, (with some results from thekinematics of finite rotations):$\frac{\partial{r\left( {p,s} \right)}}{\partial q_{j}} = {\frac{\partial{r(p)}}{\partial q_{j}} - {\left( {{{{}_{\quad}^{}{}_{\quad}^{}}(p)}{\overset{\sim}{\rho}\left( {p,s} \right)}} \right){\lambda (p)}}}$

[0255] Augmenting this equation with the additional equation λ(p,s)=λ(p)and defining w, the spatial derivative

[0256] the result is: $w\overset{\Delta}{=}\begin{bmatrix}{\lambda \quad} \\\frac{\partial r}{\partial q}\end{bmatrix}$ ${w\left( {p,s} \right)} = {\begin{bmatrix}{\lambda \left( {p,s} \right)} \\\frac{\partial{r\left( {p,s} \right)}}{\partial q_{j}}\end{bmatrix} = {{\varphi \left( {p,s} \right)}^{*}\begin{bmatrix}{\lambda (p)} \\\frac{\partial{r(p)}}{\partial q_{j}}\end{bmatrix}}}$

[0257] The vector λ can be interpreted as generating a rate of change oforientation for each body. It is a field quantity, in the sense that itcan potentially vary at each point in space. For rigid bodies undergoingpure rotations (without deformation), it is constant The second term,T₂(k), is given by:${T_{2}(k)} = {\sum\limits_{{i\varepsilon}\quad k}{\frac{\partial{\varphi \left( {k,i} \right)}}{\partial q_{j}}{T\left( {k,i} \right)}}}$

[0258] and involves the original spatial loads T(k) and derivative ofthe operator φ(k,i):

[0259] Thus ${\varphi \left( {k,i} \right)} = {\begin{bmatrix}{\underset{\underset{\_}{\_}}{E}}_{3} & {\overset{\sim}{\rho}\left( {k,i} \right)} \\{\underset{\underset{\_}{\_}}{0}}_{3} & {\underset{\underset{\_}{\_}}{E}}_{3}\end{bmatrix}\begin{bmatrix}{{{}_{\quad}^{}{}_{\quad}^{}}(k)} & {\underset{\underset{\_}{\_}}{0}}_{3} \\{\underset{\underset{\_}{\_}}{0}}_{3} & {{{}_{\quad}^{}{}_{\quad}^{}}(k)}\end{bmatrix}}^{*}$${{T_{2}(k)} = {\sum\limits_{{i\varepsilon}\quad k}{\begin{bmatrix}{{{}_{\quad}^{}{}_{\quad}^{}}(k)} & {{\overset{\sim}{\rho}\left( {k,i} \right)}{{{}_{\quad}^{}{}_{\quad}^{}}(k)}} \\{\underset{\underset{\_}{\_}}{0}}_{3} & {{{}_{\quad}^{}{}_{\quad}^{}}(k)}\end{bmatrix}^{*}{T\left( {k,i} \right)}}}}$

[0260] where

^(N) dC ^(k)(k)=^(N) dC ^(k)(i)^(i) C ^(k)(k)+^(N) C ^(k)(i)^(i) dC^(k)(k)

[0261] is computed recursively from the base body outward and${{{{}_{}^{}{}_{}^{}}(k)} = {{\frac{\partial^{i}C^{k}}{\partial{q(k)}}{{dq}(k)}\quad k} = 1}},{\ldots \quad n}$

[0262] where dq(k) is defined in the next section, and$\frac{\partial^{i}{C^{k}(k)}}{\partial q_{k}},$

[0263] the partial derivative of the interbody direction cosine matrixis a function of the joint type connecting the bodies:${\text{pin:}\quad \frac{\partial^{i}{C^{k}(k)}}{\partial q_{k}}} = {{{- {\underset{\_}{\underset{\_}{E}}}_{3}}{\sin \left( q_{k} \right)}} + {\overset{\sim}{\lambda}\quad {\cos \left( q_{k} \right)}} + {{\lambda\lambda}^{*}\quad \sin \quad \left( q_{k} \right)}}$${\text{slider:}\quad \frac{\partial^{i}{C^{k}(k)}}{\partial q_{k}}} = {\underset{\_}{\underset{\_}{0}}}_{3}$${\text{ball~~and~~free:}\quad \frac{\partial^{i}{C^{k}(k)}}{\partial ɛ_{1}}} = {2\begin{bmatrix}0 & ɛ_{2} & ɛ_{3} \\ɛ_{2} & {{- 2}ɛ_{1}} & {- ɛ_{4}} \\ɛ_{3} & ɛ_{4} & {{- 2}ɛ_{1}}\end{bmatrix}}$$\frac{\partial^{i}{C^{k}(k)}}{\partial ɛ_{2}} = {2\begin{bmatrix}{{- 2}ɛ_{2}} & ɛ_{1} & ɛ_{4} \\ɛ_{1} & 0 & ɛ_{3} \\{- ɛ_{4}} & ɛ_{3} & {{- 2}ɛ_{2}}\end{bmatrix}}$$\frac{\partial^{i}{C^{k}(k)}}{\partial ɛ_{3}} = {2\begin{bmatrix}{{- 2}ɛ_{3}} & {- ɛ_{4}} & ɛ_{1} \\ɛ_{4} & {{- 2}ɛ_{3}} & ɛ_{2} \\ɛ_{1} & ɛ_{2} & 0\end{bmatrix}}$$\frac{\partial^{i}{C^{k}(k)}}{\partial ɛ_{4}} = {2\begin{bmatrix}0 & {- ɛ_{3}} & ɛ_{2} \\ɛ_{3} & 0 & {- ɛ_{1}} \\{- ɛ_{2}} & ɛ_{1} & 0\end{bmatrix}}$

[0264] In the next Section the Force Jacobians are combined with theInertia Force Jacobians to finally form the Jacobian of the ResidualRoutine.

[0265] The Residual Jacobian

[0266] The previous Section describes the formation of the ForceJacobian $\frac{\partial{T(k)}}{\partial{w(p)}},$

[0267] which must be coupled to the spatial displacement gradient, inorder to form the derivative of the spatial forces. This sectiondescribes the formation of the spatial displacement gradient and theformation of the Jacobian of the residual routine.

[0268] The recursive algorithms for computing the entire Jacobian aredescribed. The Jacobian algorithm is not actually set up to compute theJacobian. As is typical of automatic differentiation routines, itcomputes the matrix vector product J_(uq)dq+J_(uu)du for arbitrarypassed-in values of the vectors dq and du. In practice, to compute theJacobian, the “Jacobian Routine” is effectively called repeatedly with aseries of Boolean vectors (a vector with one entry set to 1 and allother entries set to zero.) Each call generates the corresponding columnof the Jacobian. Note that some of the steps have already been or arebeing computed for the Residual Form method or the Direct Form method(the Forward Dynamics Calculations), but are reproduced here forclarity.

[0269] 1. Given (q,u) compute z

−M⁻¹ρ_(u)(q,u,0) using the Direct Form method. Also set {dot over (u)}=zand recompute A(k), then ρ_(u), which recomputes {circumflex over(T)}(k).

[0270] 2. Perform contraction to compute the fully row- andcolumn-reduced Force Jacobian, $\frac{\partial{T(k)}}{\partial{w(p)}}$

[0271] as described in section, “Interface Contraction”:$\frac{\partial T}{\partial w} = {\Phi \frac{\partial T}{\partial r}\Phi^{*}}$

[0272] Steps 3 through 10 below are used to fill the columns of J_(uq):

[0273] 3. Compute the analytic Jacobian partitions of the {dot over (q)}terms: $J_{qq} = {\frac{\partial{W(q)}}{\partial q}u}$ J_(qu) = W

[0274] using joint routines similar to those needed for the FirstKinematics Calculations.

[0275] 4. Compute q derivatives of position quantities and terms forspatial gradient:

[0276] Previously described methods were used to fill in certainjoint-specific fields. These quantities consisted of ^(i)C^(k)(k), theinterbody direction cosine matrices, r^(Q) ^(_(i)) ^(Q) ^(_(k)) (k), thespanning vector for each body, and H(k), the joint map for each body'sinboard joint. To refer to the derivative of each of these quantities, aprefix ‘d’ is added to the symbol name to make this referencegenerically. Thus, ^(i)dC^(k)(k) means the derivative of the directioncosine matrix ^(i)C^(k)(k).

[0277] Each interbody direction cosine matrix (and all thejoint-specific) quantities depend only on the generalized coordinates ofan individual joint. Thus, ^(i)dC^(k)(k) is nonzero only when thederivative is taken with respect to any of the coordinates for body k.To properly ‘seed’ the recursions being studying, a vector dq is passedin to the routine. For Jacobian computation we set one entry is set to1, and all the other entries to 0. Then, the needed preliminaryquantities are generated by a typical loop:${{{{}_{}^{}{}_{}^{}}(k)} = {{\frac{\partial^{i}C^{k}}{\partial{q(k)}}{{dq}(k)}\quad k} = 1}},{\ldots \quad n}$

[0278] The partial derivatives of the direction cosine matrices aregenerated analytically and displayed in the section, “InterfaceContraction” above. These partial derivatives do not depend upon theparticular column of the Jacobian that is being computed. Setting aparticular entry of dq to 1, and all the rest to 0 has the effect ofannihilating the correct subset of the seed quantities.$\frac{\partial{r^{Q_{i}Q_{k}}(k)}}{\partial q_{k}},$

[0279] the partial derivative of the interbody spanning vector is givenby${\frac{\partial{r^{Q_{i}Q_{k}}(k)}}{\partial q_{k}} = {\lambda (k)}},\text{slider}$${\frac{\partial{r^{Q_{i}Q_{k}}(k)}}{\partial q_{k}} = {\underset{\_}{\underset{\_}{0}}}_{3{x4}}},\text{ball}$${\frac{\partial{r^{Q_{i}Q_{k}}(k)}}{\partial q_{k}} = {\underset{\_}{0}}_{3}},\text{pin}$${\frac{\partial{r^{Q_{i}Q_{k}}(k)}}{\partial q_{k}} = \begin{bmatrix}0 & 0 & 0 & 0 & 1 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 1 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 1\end{bmatrix}},\text{free}$

[0280] λ(k) here refers to the body's sliding axis which connects it toits parent body. $\frac{\partial{H(k)}}{\partial q_{k}},$

[0281] the partial derivative of the joint map is${\frac{\partial{H(k)}}{\partial q_{k}} = {\underset{\_}{0}}_{6}^{*}},{pin},{slider}$${\frac{\partial{H(k)}}{\partial q_{k}} = \begin{bmatrix}{\underset{\underset{\_}{\_}}{0}}_{3} & {\underset{\underset{\_}{\_}}{0}}_{3}\end{bmatrix}},{ball}$${\frac{\partial{H(k)}}{\partial q_{k}} = \begin{bmatrix}{\underset{\underset{\_}{\_}}{0}}_{3} & {\underset{\underset{\_}{\_}}{0}}_{3} \\{\underset{\underset{\_}{\_}}{0}}_{3} & \frac{\partial^{i}{C_{\quad}^{k}(k)}}{\partial q_{k}}\end{bmatrix}},{free}$

[0282] With the above definitions of the partial derivatives, therecursions are seeded with the following loop:for  k = 1  to  nbod${{{}_{}^{}{}_{\quad}^{}}(k)} = {\frac{\partial^{i}{C_{\quad}^{k}(k)}}{\partial q_{k}}{{dq}(k)}}$${{dr}^{Q_{i}Q_{k}}(k)} = {\frac{\partial{r^{Q_{i}Q_{k}}(k)}}{\partial q_{k}}{{dq}(k)}}$${{dH}(k)} = {\frac{\partial{H(k)}}{\partial q_{k}}{{dq}(k)}}$ end

[0283] After execution of these loops, all bodies have ^(i)dC^(k)(k),dr^(Q) ^(_(i)) ^(Q) ^(_(k)) (k), and dH(k), their interbody derivativequantities available.

[0284] One new quantity needed in the spatial displacement gradientcomputation is also computed. This is λ(k) from the section on InterfaceContraction, the rotation axis that generates the rate of change oforientation for each body outboard of a perturbed joint. Here, thisvariable is given the symbol dθ(k), the differential rotation axis foreach body is

[0285] for k=1 to nbod

[0286] dθ(k)=λ(k) dq(k), pin

[0287] dθ(k)= 0 ₃, slider

[0288] dθ(k)=not needed for ball, free

[0289] end

[0290] Since arbitrary perturbations to a set of Euler parameters do notproduce a pure rotation, column contraction cannot be used whencomputing the corresponding column of the Jacobian. The row-reducedForce Jacobian can still (and must) be used.

[0291] After seeding the recursions, ^(N)dC^(k)(k), dr^(OQ) ^(_(k)) (k),^(i)dφ^(k)(k), dλ(k) are computed:   (1) =   (1)dr^(OQ_(k))(1) = dr^(Q_(i)Q_(k))(1)${{{}_{}^{}{}_{\quad}^{}}(1)} = {\underset{\underset{\_}{\_}}{0}}_{6}$dλ(1) = dθ(1) for  k = 2  to  nbod  (k) =   (i)  (k) +   (i)  (k)dr^(OQ_(k))(k) = dr^(OQ_(k))(i) +   (i)r^(Q_(i)Q_(k))(k) +   (i)dr^(Q_(i)Q_(k))(k)${{{}_{}^{}{}_{\quad}^{}}(k)} = \begin{bmatrix}a & b \\c & d\end{bmatrix}$ dλ(k) =   ^(k*)(k)dλ(i) + dθ(k) end where${a = {{{}_{}^{}{}_{\quad}^{}}(k)}},{b = {{d{{\overset{\sim}{r}}^{Q_{i}Q_{k}}(k)}{{{}_{}^{}{}_{\quad}^{}}(k)}} + {{{}_{}^{}{}_{\quad}^{}}(k)}}},{c = {\underset{\underset{\_}{\_}}{0}}_{3}},{d = {{{}_{}^{}{}_{\quad}^{}}(k)}}$

[0292] 5. Compute q derivatives of velocity:

[0293] A loop that computes the rate of change of joint velocity due tochange in joint angle starts the process:

[0294] for k=1 to nbod

[0295]^(i)dν^(k)(k)=dH*(k)u(k)

[0296] end

[0297] This quantity is the rate of change of joint velocity due tochange in joint angle. Obviously, it is nonzero only for joints whosemap contains coordinate dependence. For free joints, the generalizedspeeds produce relative linear velocity that depends upon the jointorientation.

[0298] After computing ^(i)dν^(k)(k), dV(k), the derivative of thespatial velocity of each body, is computed. This is done by thefollowing loop:

[0299] dV(1)=^(i)dν^(k)(1)

[0300] for k=2 to nbod

[0301] dV(k)=^(i)dφ^(k)*V(i)+^(i)φ^(k*)dV(i)+^(i)dν^(k)(k)

[0302] end

[0303] 6. Couple Force Jacobian to spatial displacement gradient tocompute T₁(k) for  k = 1  to  nbod${T_{1}(k)} = {\sum\limits_{p = 1}^{nbod}{\frac{\partial{T(k)}}{\partial{w(p)}}\frac{\partial{w(p)}}{\partial q_{j}}}}$

[0304] end

[0305] 7. Compute second term of the Force Jacobian T₂(k) and append toT₁(k):

[0306] for k=1 to nbod${{dT}(k)} = {{T_{1}(k)} + {\sum\limits_{{i\varepsilon}\quad k}{\begin{bmatrix}{{{}_{}^{}{}_{\quad}^{}}(k)} & {{\overset{\sim}{\rho}\left( {k,i} \right)}^{N}{{dC}^{k}(k)}} \\{\underset{\underset{\_}{\_}}{0}}_{3} & {{{}_{}^{}{}_{\quad}^{}}(k)}\end{bmatrix}^{*}{T\left( {k,i} \right)}}}}$

[0307] end

[0308] 8. Compute q derivatives of acceleration-related terms:

[0309] Again the process starts with a loop that computes^(i)da^(k)(k)=dH*(k){dot over (u)}(k):

[0310] for k=1 to nbod

[0311]^(i)da^(k)(k)=dH*(k){dot over (u)}(k)

[0312] end

[0313] This is the change in joint acceleration due to a change incoordinate. Then, dA(k), the derivative of the spatial acceleration ofeach body is computed. dA(1) =   (1)$\left. {{{}_{}^{}{}_{\quad}^{}}(k)}\rightarrow\begin{bmatrix}{{{}_{}^{}{}_{\quad}^{}}(k)} \\{{{}_{}^{}{}_{\quad}^{Qk}}(k)}\end{bmatrix} \right.,\left. {{{}_{}^{}{}_{\quad}^{}}(k)}\rightarrow\begin{bmatrix}{{{}_{}^{}{}_{\quad}^{}}(k)} \\{{{}_{}^{}{}_{\quad}^{Qk}}(k)}\end{bmatrix} \right.$ $\left. {V(k)}\rightarrow\begin{bmatrix}{{{}_{}^{}{}_{\quad}^{}}(k)} \\{{{}_{}^{}{}_{\quad}^{Qk}}(k)}\end{bmatrix} \right.,\left. {{dV}_{\quad}(k)}\rightarrow\begin{bmatrix}{{{}_{}^{}{}_{\quad}^{}}(k)} \\{{{}_{}^{}{}_{\quad}^{Qk}}(k)}\end{bmatrix} \right.$

[0314] where → means the 6 vector is decomposed into two 3 vectors, fork=2 to nbod${dt1}\overset{\Delta}{=}{{{{{}_{}^{}{}_{\quad}^{k*}}(k)}\left( {{{{}_{}^{}{}_{}^{}}(i)} \times \left( {{{{}_{}^{}{}_{}^{}}(i)} \times {r^{Q_{i}Q_{k}}(k)}} \right)} \right)} + {{{{}_{}^{}{}_{}^{k*}}(k)}\begin{pmatrix}{\left( {{{{}_{}^{}{}_{}^{}}(i)} \times \left( {{{{}_{}^{}{}_{}^{}}(i)} \times {r^{Q_{i}Q_{k}}(k)}} \right)} \right) +} \\{\left( {{{{}_{}^{}{}_{}^{}}(i)} \times \left( {{{{}_{}^{}{}_{}^{}}(i)} \times {r^{Q_{i}Q_{k}}(k)}} \right)} \right) +} \\\left( {{{{}_{}^{}{}_{}^{}}(i)} \times \left( {{{{}_{}^{}{}_{}^{}}(i)} \times {{dr}^{Q_{i}Q_{k}}(k)}} \right)} \right)\end{pmatrix}}}$${da}\overset{\Delta}{=}{{{{{}_{}^{}{}_{}^{k*}}(k)}{A(i)}} + {{{{}_{}^{}{}_{}^{k*}}(k)}{{dA}(i)}} + \begin{bmatrix}{\underset{\underset{\_}{\_}}{0}}_{3} \\{dt1}\end{bmatrix}}$${\omega \quad j}\overset{\Delta}{=}{{{{}_{}^{}{}_{}^{k*}}(k)}{{{}_{}^{}{}_{}^{}}(i)}}$${{d\omega}\quad j}\overset{\Delta}{=}{{{{{}_{}^{}{}_{}^{k*}}(k)}{{{}_{}^{}{}_{}^{}}(i)}} + {{{{}_{}^{}{}_{}^{k*}}(k)}{{{}_{}^{}{}_{}^{}}(i)}}}$${dt2}\overset{\Delta}{=}{{{d\omega}\quad j \times {{{}_{}^{}{}_{}^{}}(k)}} + {\omega \quad j \times {{{}_{}^{}{}_{}^{}}(k)}}}$${dt3}\overset{\Delta}{=}{{2{d\omega}\quad j \times {{{}_{}^{}{}_{}^{Qk}}(k)}} + {2\omega \quad j \times {{{}_{}^{}{}_{}^{Qk}}(k)}}}$${{dA}(k)} = {{da} + \begin{bmatrix}{dt2} \\{dt3}\end{bmatrix} + {{{}_{}^{}{}_{}^{}}(k)}}$ end

[0315] The symbols introduced here with

are meant to be temporary variables not needed after computation ofdA(k).

[0316] After computing the spatial acceleration derivatives, thecomputation of d{circumflex over (T)}(k), the spatial inertia forcederivatives, is performed:

[0317] for k=1 to nbod

[0318] dν1

^(N)dω^(k)(k)×I _(Q) _(k) (k)·^(N)ω^(k)+^(N)ω^(k)(k)×I _(Q) _(k)(k)·^(N)dω^(k)

[0319] dν2

^(N)dω^(k)(k)×(^(N)ω^(k)(k)×p(k))+^(N)ω^(k)(k)×(^(N)dω^(k)(k)×p(k))${d{\hat{T}(k)}} = {{{M(k)}{{dA}(k)}} + \begin{bmatrix}{dv1} \\{dv2}\end{bmatrix} - {{dT}(k)}}$

[0320] end

[0321] 9. Compute dρ(k), the joint residual derivative for body k:

[0322] for k=nbod to 1

[0323] dρ(k)=dH(k){circumflex over (T)}(k)+H(k)d{circumflex over(T)}(k)−dσ(k)

[0324] i=inb(k)

[0325] if i>0

[0326] d{circumflex over (T)}(i)=d{circumflex over(T)}(i)+^(i)dφ^(k)(k){circumflex over (T)}(k)+^(i)φ^(k)(k)d{circumflexover (T)}(k)

[0327] end

[0328] end

[0329] After executing this routine, the values stored in dρ(k) are thenew column of the Residual Jacobian$\frac{\partial\quad}{\partial q}{{\rho_{u}\left( {q,u,z} \right)}.}$

[0330] 10. Back-solve the $\frac{\partial\rho}{\partial q}$

[0331] result of previous step with the mass-matrix to obtain thedesired $\frac{\partial\overset{.}{u}}{\partial q}:$

$\frac{\partial\overset{.}{u}}{\partial q} = \quad {{- M^{- 1}}\frac{\partial\rho}{\partial q}}$

[0332] The back-solve operation is accomplished in the Direct Formmethod routine by processing a residual vector into a {dot over (u)}vector. The Second Kinematics Calculations only needs to be performedonce for the whole Jacobian, since the back-solves are done at thenominal value of the state. In fact, the Second Kinematics routine musthave been called in Step 1 while computing z, so the variables shouldstill be cached.

[0333] Steps 11 through 13 below are used to fill the columns of J_(uu):

[0334] 11. Compute u derivatives of velocity:

[0335] This routine takes a passed-in vector du and computes^(i)dν^(k)(k)=H*(k)du(k). Then, dV(k), the derivative of the spatialvelocity of each body, is computed:

[0336] dV(1)=^(i)dν^(k)(1)

[0337] for k=2 to nbod

[0338] dV(k)=^(i)φ^(k)*(k)dV(i)+^(i)dν^(k)(k)

[0339] end

[0340] 12. Compute the velocity-induced derivative d{circumflex over(T)}(k). As presented here, the routine is specialized for the case ofno velocity dependent external loads. The surviving terms are those dueto changes in inertia forces alone. Even if there were changes inexternal loads, it would only be required to include the contribution ofdT(k) as before. ${{dA}(1)} = \begin{bmatrix}{\underset{\underset{\_}{\_}}{0}}_{3} \\{\underset{\underset{\_}{\_}}{0}}_{3}\end{bmatrix}$ for  k = 2  to  nbod${dt1}\overset{\Delta}{=}{{{{}_{}^{}{}_{}^{k*}}(k)}\begin{pmatrix}{\left( {{{{}_{}^{}{}_{}^{}}(i)} \times \left( {{{{}_{}^{}{}_{}^{}}(i)} \times {r^{Q_{i}Q_{k}}(k)}} \right)} \right) +} \\\left( {{{{}_{}^{}{}_{}^{}}(i)} \times \left( {{{{}_{}^{}{}_{}^{}}(i)} \times {r^{Q_{i}Q_{k}}(k)}} \right)} \right)\end{pmatrix}}$${da}\overset{\Delta}{=}{{{{{}_{}^{}{}_{}^{k*}}(k)}{{dA}(i)}} + \begin{bmatrix}{\underset{\_}{0}}_{3} \\{dt1}\end{bmatrix}}$${\omega \quad j}\overset{\Delta}{=}{{{{}_{}^{}{}_{}^{k*}}(k)}{{{}_{}^{}{}_{}^{}}(i)}}$${{d\omega}\quad j}\overset{\Delta}{=}{{{{}_{}^{}{}_{}^{k*}}(k)}{{{}_{}^{}{}_{}^{}}(i)}}$${dt2}\overset{\Delta}{=}{{{d\omega}\quad j \times {{{}_{}^{}{}_{}^{}}(k)}} + {\omega \quad j \times {{{}_{}^{}{}_{}^{}}(k)}}}$${dt3}\overset{\Delta}{=}{{2{d\omega}\quad j \times {{{}_{}^{}{}_{}^{Qk}}(k)}} + {2\omega \quad j \times {{{}_{}^{}{}_{}^{Qk}}(k)}}}$${{dA}(k)} = {{da} + \begin{bmatrix}{dt2} \\{dt3}\end{bmatrix}}$

[0341] After computing the spatial acceleration derivatives,d{circumflex over (T)}(k), the spatial inertia force derivatives, iscomputed: for  k = 1  to  nbod${dv1}\overset{\Delta}{=}{{{\,^{N}d}\quad {\omega^{k}(k)} \times {{{\underset{\underset{\_}{\_}}{I}}_{Q_{k}}(k)} \cdot {{}_{}^{}{}_{}^{}}}} + {{{{}_{}^{}{}_{}^{}}(k)} \times {{{\underset{\underset{\_}{\_}}{I}}_{Q_{k}}(k)} \cdot {\,^{N}d}}\quad \omega^{k}}}$${dv2}\overset{\Delta}{=}{{{\,^{N}d}\quad {\omega^{k}(k)} \times \left( {{{{}_{}^{}{}_{}^{}}(k)} \times {p(k)}} \right)} + {{{{}_{}^{}{}_{}^{}}(k)} \times \left( {{\,^{N}d}\quad {\omega^{k}(k)} \times {p(k)}} \right)}}$${d{\hat{T}(k)}} = {{{M(k)}{{dA}(k)}} + \begin{bmatrix}{dv1} \\{dv2}\end{bmatrix}}$ end

[0342] 13. Compute dρ(k), the joint residual derivative for body k:

[0343] for k=nbod to 1

[0344] dρ(k)=H(k)d{circumflex over (T)}(k)−dσ(k)

[0345] i=inb(k)

[0346] if i>0

[0347] d{circumflex over (T)}(i)=d{circumflex over(T)}(i)+^(i)φ^(k)(k)d{circumflex over (T)}(k)

[0348] end

[0349] end

[0350] After executing this routine the values stored in dρ(k) are thenew column of the Residual Jacobian$\frac{\partial\quad}{\partial u}{{\rho_{u}\left( {q,u,z} \right)}.}$

[0351] 14. Back-solve the $\frac{\partial\rho}{\partial u}$

[0352] result of previous step with the mass-matrix to obtain thedesired $\frac{\partial\overset{.}{u}}{\partial u}:$

$\frac{\partial\overset{.}{u}}{\partial u} = \quad {{- M^{- 1}}\frac{\partial\rho}{\partial u}}$

[0353] The back-solve operation is accomplished in the Direct Methodroutine by processing a residual vector into a {dot over (u)} vector.The Second Kinematics Calculations only need to be performed once, sincethe back-solves are done at the nominal value of the state. In fact, theSecond Kinematics routine must have been called in Step 1 whilecomputing z, so the variables should still be cached.

[0354] The above steps complete the computation of the analytic Jacobianas long as the forces only have dependence on q. This accommodates theclassical situation where all atomic forces are derivable from apotential function. To accommodate velocity-dependent forces, such assimple viscous damping, some of the above steps need to be modified asfollows:

[0355] In Step 2 above, we also need to compute the contracted velocityJacobian $\frac{\partial{T(k)}}{\partial{\overset{.}{r}(k)}},$

[0356] which is block diagonal, must also be computed.

[0357] In Step 6 above, the computation of T₁(k) must be augmented withthe contracted velocity Jacobian:

[0358] for k=1 to nbod${T_{1}(k)} = {{\sum\limits_{p = 1}^{nbod}{\frac{\partial{T(k)}}{\partial{w(p)}}\frac{\partial{w(p)}}{\partial q_{j}}}} + {\sum\limits_{{i\varepsilon}\quad k}{\frac{\partial{T(k)}}{\partial{\overset{.}{r}\left( {k,i} \right)}}\frac{\partial{\overset{.}{r}\left( {k,i} \right)}}{\partial q_{j}}}}}$

[0359] end

[0360] where$\frac{\partial{\overset{.}{r}\left( {k,i} \right)}}{\partial q_{j}} = {{{{{}_{}^{}{}_{}^{}}(k)}\left\lbrack {{{{}_{}^{}{}_{}^{Qk}}(k)} + {{{{}_{}^{}{}_{}^{}}(k)} \times {\rho \left( {k,i} \right)}}} \right\rbrack} + {{{{}_{}^{}{}_{}^{}}(k)}\left\lbrack {{{{}_{}^{}{}_{}^{Qk}}(k)} + {{{{}_{}^{}{}_{}^{}}(k)} \times {\rho \left( {k,i} \right)}}} \right\rbrack}}$

[0361] A step is then added after Step 11, which is called Step 11a.This new step computes dT(k):

[0362] where${{dT}(k)} = {\sum\limits_{i \in k}{\frac{\partial{T(k)}}{\partial{\overset{.}{r}\left( {k,i} \right)}}\frac{\partial{\overset{.}{r}\left( {k,i} \right)}}{\partial u_{j}}}}$where$\frac{\partial{\overset{.}{r}\left( {k,i} \right)}}{\partial u_{j}} = {{{{}_{}^{}{}_{}^{}}(k)}\left\lbrack {{{{}_{}^{}{}_{}^{Qk}}(k)} + {{\,^{N}d}\quad {\omega^{k}(k)} \times {\rho \left( {k,i} \right)}}} \right\rbrack}$

[0363] While executing Step 12 above, the last loop for d{circumflexover (T)}(k) is modified by subtracting the velocity-dependent forcederivative dT(k): for  k = 1  to  nbod${dv1}\overset{\Delta}{=}{{{\,^{N}d}\quad {\omega^{k}(k)} \times {{{\underset{\_}{\underset{\_}{I}}}_{Q_{k}}(k)} \cdot {{}_{}^{}{}_{}^{}}}} + {{{{}_{}^{}{}_{}^{}}(k)} \times {{{\underset{\_}{\underset{\_}{I}}}_{Q_{k}}(k)} \cdot {\,^{N}d}}\quad \omega^{k}}}$${dv2}\overset{\Delta}{=}{{{\,^{N}d}\quad {\omega^{k}(k)} \times \left( {{{{}_{}^{}{}_{}^{}}(k)} \times {p(k)}} \right)} + {{{{}_{}^{}{}_{}^{}}(k)} \times \left( {{\,^{N}d}\quad {\omega^{k}(k)} \times {p(k)}} \right)}}$${d{\hat{T}(k)}} = {{{M(k)}{{dA}(k)}} + \begin{bmatrix}{dv1} \\{dv2}\end{bmatrix} - {{dT}(k)}}$

[0364] end

[0365] The rest of the Steps remain the same.

[0366]FIG. 6 summarizes the operational steps of the Analytic Jacobianmethod, which has been described in detail above.

[0367]FIG. 7 shows a plot of the accuracy of the numerical Jacobianversus the accuracy of the analytic Jacobian for an exemplary MD system.In the best case in which the perturbation was perfectly selected, thedigits of accuracy for the generalized coordinates (q) and generalizedspeeds (u) from the numerical Jacobian, illustrated by line 152, werestill only half that of the analytic Jacobian, illustrated by line 150.

[0368] Additional Embodiments

[0369] The present invention has many embodiments besides the examplesdescribed above. The list below has other embodiments and applications:

[0370] Order of Forces Included in Jacobian

[0371] Any order of the forces to be included in the Jacobian, include,but not limited to Order(N), Order(N²), Order(N³), and Order(N⁴). Anexample of an Order(N) force field would be an electrostatic force fieldusing fast multi-pole expansion methods (see, for example, Greengaard,The Rapid Evaluation of Potential Fields in Particle Systems,Massachusetts Institute of Technology Dissertation, 1988) rather thandirect method which is Order(N²).

[0372] Analytic Jacobian for Direct Form

[0373] When the governing equations are in Direct Form, the so-called“forward dynamics” form of the equations is obtained. In this form, theequations process a state vector and applied efforts and generate theacceleration at each of the joints modeled in the system.

{dot over (u)}=M ⁻¹(ƒ)

[0374] The Jacobian then represents the partial derivatives of theaccelerations with respect to elements of the state vector. Thepreferred embodiment shows several algorithmic methods for computationof these partial derivatives. The methods are exact and do not utilizenumerical approximations to form derivatives.

[0375] The Direct Form produces the {dot over (u)} partitions of theJacobian:

[0376] and$J_{uq} = \frac{\partial\left( {M^{- 1}(f)} \right)}{\partial q}$ and$J_{uu} = \frac{\partial\left( {M^{- 1}(f)} \right)}{\partial u}$

[0377] by using an algorithmic counterpart to the function whichcomputes the {dot over (u)} function.

What is claimed is:
 1. A method of modeling the behavior of a molecule,comprising selecting a torsion angle, rigid multibody model for saidmolecule, said model having equations of motion; selecting an implicitintegrator; and generating an analytic Jacobian for said implicitintegrator to integrate said equations of motion so as to obtaincalculations of said behavior of said molecule.
 2. The method of claim 1wherein said analytic Jacobian is derived from an analytic Jacobian ofthe Residual Form of the equations of motion.
 3. The method of claim 2wherein said analytic Jacobian J comprises ${J = {\begin{pmatrix}\frac{\partial\overset{.}{q}}{\partial q} & \frac{\partial\overset{.}{q}}{\partial u} \\\frac{\partial\overset{.}{u}}{\partial q} & \frac{\partial\overset{.}{u}}{\partial u}\end{pmatrix}\overset{\Delta}{=}\begin{pmatrix}J_{qq} & J_{qu} \\J_{uq} & J_{uu}\end{pmatrix}}};\text{and}$$J_{qq} = {\frac{\partial\overset{.}{q}}{\partial q} = {{\frac{\partial\left( {W\quad u} \right)}{\partial q}\quad \text{and}\quad J_{qu}} = {\frac{\partial\overset{.}{q}}{\partial u} = W}}}$$J_{uq} = {\frac{\partial\overset{.}{u}}{\partial q} = {{{- M^{- 1}}\frac{\partial{\rho_{u}\left( {q,u,z} \right)}}{\partial q}\quad \text{and}\quad J_{uu}} = {\frac{\partial\overset{.}{u}}{\partial u} = {{- M^{- 1}}\frac{\partial{\rho_{u}\left( {q,u,z} \right)}}{\partial u}}}}}$

where q are the generalized coordinates, u are the generalized speeds, Wis a joint map matrix and M is the mass matrix and ρ_(u) is the dynamicresidual of the equations of motion, and z is −M⁻¹ρ_(u)(q,u,0).
 4. Themethod of claim 3 wherein said implicit integrator selecting stepcomprises an L-stable integrator.
 5. A method of simulating the behaviorof a physical system, comprising modeling said physical system with atorsion angle, rigid multibody model, said model having equations ofmotion; and integrating said equations of motion with an implicitintegrator; said implicit integrator having an analytic Jacobian toobtain calculations of said behavior of said physical system.
 6. Themethod of claim 5 wherein said analytic Jacobian is derived from ananalytic Jacobian of the Residual Form of the equations of motion. 7.The method of claim 6 wherein said analytic Jacobian J comprises${J = {\begin{pmatrix}\frac{\partial\overset{.}{q}}{\partial q} & \frac{\partial\overset{.}{q}}{\partial u} \\\frac{\partial\overset{.}{u}}{\partial q} & \frac{\partial\overset{.}{u}}{\partial u}\end{pmatrix}\overset{\Delta}{=}\begin{pmatrix}J_{qq} & J_{qu} \\J_{uq} & J_{uu}\end{pmatrix}}};\text{and}$$J_{qq} = {\frac{\partial\overset{.}{q}}{\partial q} = {{\frac{\partial\left( {W\quad u} \right)}{\partial q}\quad \text{and}\quad J_{qu}} = {\frac{\partial\overset{.}{q}}{\partial u} = W}}}$$J_{uq} = {\frac{\partial\overset{.}{u}}{\partial q} = {{{- M^{- 1}}\frac{\partial{\rho_{u}\left( {q,u,z} \right)}}{\partial q}\quad \text{and}\quad J_{uu}} = {\frac{\partial\overset{.}{u}}{\partial u} = {{- M^{- 1}}\frac{\partial{\rho_{u}\left( {q,u,z} \right)}}{\partial u}}}}}$

where q are the generalized coordinates, u are the generalized speed, Wis a joint map matrix and M is the mass matrix and ρ_(u) is the dynamicresidual of the equations of motion, and z is −M⁻¹ρ_(u)(q,u,0).
 8. Themethod of claim 7 wherein said implicit integrator comprises an L-stableintegrator.
 9. Computer code for simulating the behavior of a molecule,said code comprising a first module for a torsion angle, rigid multibodymodel of said molecule, said model having equations of motion; and asecond module for an implicit integrator to integrate said equations ofmotion with an analytic Jacobian to obtain calculations of said behaviorof said molecule.
 10. The computer code of claim 9 wherein said analyticJacobian is derived from an analytic Jacobian of the Residual Form ofthe equations of motion.
 11. The computer code of claim 10 wherein saidanalytic Jacobian J comprises ${J = {\begin{pmatrix}\frac{\partial\overset{.}{q}}{\partial q} & \frac{\partial\overset{.}{q}}{\partial u} \\\frac{\partial\overset{.}{u}}{\partial q} & \frac{\partial\overset{.}{u}}{\partial u}\end{pmatrix}\overset{\Delta}{=}\begin{pmatrix}J_{qq} & J_{qu} \\J_{uq} & J_{uu}\end{pmatrix}}};\text{and}$$J_{qq} = {\frac{\partial\overset{.}{q}}{\partial q} = {{\frac{\partial\left( {W\quad u} \right)}{\partial q}\quad \text{and}\quad J_{qu}} = {\frac{\partial\overset{.}{q}}{\partial u} = W}}}$$J_{uq} = {\frac{\partial\overset{.}{u}}{\partial q} = {{{- M^{- 1}}\frac{\partial{\rho_{u}\left( {q,u,z} \right)}}{\partial q}\quad \text{and}\quad J_{uu}} = {\frac{\partial\overset{.}{u}}{\partial u} = {{- M^{- 1}}\frac{\partial{\rho_{u}\left( {q,u,z} \right)}}{\partial u}}}}}$

where q are the generalized coordinates, u are the generalized speed, Wis a joint map matrix and M is the mass matrix and ρ_(u) is the dynamicresidual of the equations of motion, and z is −M⁻¹ρ_(u)(q,u,0).
 12. Thecomputer code of claim 11 wherein said implicit integrator comprises anL-stable integrator.
 13. Computer code for simulating the behavior of aphysical system, said code comprising a first module for a torsionangle, rigid multibody model of said system, said model having equationsof motion; and a second module for an implicit integrator to integratesaid equations of motion with an analytic Jacobian to obtaincalculations of said behavior of said system.
 14. The computer code ofclaim 13 wherein said analytic Jacobian is derived from an analyticJacobian of the Residual Form of the equations of motion.
 15. Thecomputer code of claim 14 wherein said analytic Jacobian J comprises${J = {\begin{pmatrix}\frac{\partial\overset{.}{q}}{\partial q} & \frac{\partial\overset{.}{q}}{\partial u} \\\frac{\partial\overset{.}{u}}{\partial q} & \frac{\partial\overset{.}{u}}{\partial u}\end{pmatrix}\overset{\Delta}{=}\begin{pmatrix}J_{qq} & J_{qu} \\J_{uq} & J_{uu}\end{pmatrix}}};\text{and}$$J_{qq} = {\frac{\partial\overset{.}{q}}{\partial q} = {{\frac{\partial\left( {W\quad u} \right)}{\partial q}\quad \text{and}\quad J_{qu}} = {\frac{\partial\overset{.}{q}}{\partial u} = W}}}$$J_{uq} = {\frac{\partial\overset{.}{u}}{\partial q} = {{{- M^{- 1}}\frac{\partial{\rho_{u}\left( {q,u,z} \right)}}{\partial q}\quad \text{and}\quad J_{uu}} = {\frac{\partial\overset{.}{u}}{\partial u} = {{- M^{- 1}}\frac{\partial{\rho_{u}\left( {q,u,z} \right)}}{\partial u}}}}}$

where q are the generalized coordinates, u are the generalized speed, Wis a joint map matrix and M is the mass matrix and ρ_(u) is the dynamicresidual of the equations of motion, and z is −M⁻¹ρ_(u)(q,u,0).
 16. Thecomputer code of claim 15 wherein said implicit integrator comprises anL-stable integrator.